We  describe  a  new  class  of  algorithms  for  the  solution  of  parabolic  partial  differential  equa¬ 
tions  (PDEs).  This  class  of  schemes  is  based  on  three  principal  observations.  First,  the 
spatial  discretization  of  parabolic  PDEs  results  in  stiff  systems  of  ordinary  differential  equa¬ 
tions  (ODEs)  in  time,  and  therefore,  requires  an  implicit  method  for  its  solution.  Spectral 
Deferred  Correction  (SDC)  methods  use  repeated  iterations  of  a  low-order  method  (e.g.  im¬ 
plicit  Euler  method)  to  generate  a  high-order  scheme.  As  a  result,  SDC  methods  of  arbitrary 
order  can  be  constructed  with  the  desired  stability  properties  necessary  for  the  solution  of 
stiff  differential  equations.  Furthermore,  for  large-scale  systems,  SDC  methods  are  more 
computationally  efficient  than  implicit  Runge-Kutta  schemes.  Second,  implicit  methods  for 
the  solution  of  a  system  of  linear  ODEs  yield  linear  systems  that  must  be  solved  on  each  it¬ 
eration.  It  is  well  known  that  the  linear  systems  constructed  from  the  spatial  discretization 
of  parabolic  PDEs  are  sparse.  In  R1,  these  linear  systems  can  be  solved  in  0(n )  operations 
where  n  is  the  number  of  spatial  discretization  nodes.  However,  in  R2,  the  straightfor¬ 
ward  spatial  discretization  leads  to  matrices  with  dimensionality  n 2  x  n2  and  bandwidth 
n.  While  fast  inversion  schemes  of  0(n3)  exist,  we  use  alternating  direction  implicit  (ADI) 
methods  to  replace  the  single  two-dimensional  implicit  step  with  two  sub-steps  where  only 
one  direction  is  treated  implicitly.  This  approach  results  in  schemes  with  computational 
cost  0(n2).  Likewise,  ADI  methods  in  R3  have  computational  cost  0(n3).  While  popular 
ADI  methods  are  low-order,  we  combine  the  SDC  methods  with  an  ADI  method  to  generate 
computationally  efficient,  high-order  schemes  for  the  solution  of  parabolic  PDEs  in  R2  and 
R3.  Third,  traditional  pseudospectral  schemes  for  the  representation  of  the  spatial  operator 
in  parabolic  PDEs  yield  differentiation  operators  with  eigenvalues  that  can  be  excessively 
large.  We  improve  on  the  traditional  approach  by  subdividing  the  entire  spatial  domain, 
constructing  bases  on  each  subdomain,  and  combining  the  obtained  discretization  with  the 
implicit  SDC  schemes.  The  resulting  schemes  are  high-order  in  both  time  and  space  and 
have  computational  cost  0{N  ■  M)  where  N  is  the  number  of  spatial  discretization  nodes 
and  M  is  the  number  of  temporal  nodes.  We  illustrate  the  behavior  of  these  schemes  with 
several  numerical  examples. 
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1  Introduction 


The  diffusion  equation  with  variable  coefficients  is  an  extremely  well-studied  subject.  It  is 
usually  written  as 


dt 


(a(x)V(f>) , 


(1.1) 


where  a  is  the  diffusion  coefficient  and  x  E  for  d  =  1,2,  3.  This  equation  arises  in  various 
physical  phenomena,  including  heat  transfer  and  fluid  dynamics,  biological  and  chemical 
reaction  processes,  financial  mathematics,  and  many  other  application  areas.  Furthermore, 
the  diffusion  equation  is  an  important  component  of  many  non-linear  partial  differential 
equations  (PDEs),  including  the  Navier-Stokes  equations.  Thus,  there  is  continued  interest 
in  constructing  accurate  and  fast  methods  for  its  solution.  Traditional  approaches  for 
computing  its  numerical  solution  include  finite  difference,  finite  element,  integral  equation, 
spectral,  and  pseudospectral  methods  (see,  e.g.,  13  0  Q21  da  n  DU  E2  m  [Ml  ED] )  • 

Finite  difference  and  finite  element  methods  tend  to  lead  to  low-order  schemes  for  the 
solution  of  (1.1)  and,  hence,  tend  to  require  an  excessive  number  of  nodes  in  both  the  tem¬ 
poral  and  spatial  discretization.  Furthermore,  high-order  finite  difference  and  finite  element 
methods  have  difficulty  maintaining  the  high-order  near  the  boundary.  As  a  consequence, 
there  is  continued  interest  in  developing  high-order  methods  in  both  time  and  space. 

For  the  case  of  constant  coefficients,  one  can  convert  (1.1)  to  an  integral  equation  and 
evaluate  the  heat  potential  (see,  e.g.,  HU  ED  [22]).  For  variable  coefficients,  some  of  the 
issues  of  spatial  discretization  mentioned  above  are  resolved  in  discontinuous  Galerkin  or 
discontinuous  spectral  element  methods  (see,  e.g.,  1231  El  [37]). 

Spatial  discretization  of  (1.1)  results  in  a  stiff  system  of  ordinary  differential  equations 
(ODEs)  in  time  and,  therefore,  requires  an  implicit  method  for  its  solution.  While  there  are 
no  implicit  A-stable  multistep  methods  of  order  greater  than  two,  there  is  a  family  of  implicit 
multistep  methods  which  are  A(a)-stable  called  backward  differentiation  formulae  (BDF). 
However,  since  there  are  no  BDFs  of  order  greater  than  six,  implicit  Runge-Kutta  schemes 
are  frequently  used  for  the  solution  of  stiff  ODEs.  On  the  other  hand,  implicit  Runge-Kutta 
schemes  are  computationally  expensive  for  large  systems.  An  alternative  approach,  which 
we  use  in  this  dissertation,  is  the  so-called  Spectral  Deferred  Correction  (SDC)  methods 
m-  To  avoid  high  cost,  SDC  methods  use  repeated  iterations  of  a  low-order  method  (e.g., 
implicit  Euler  method)  to  generate  a  high-order  scheme.  As  a  result,  SDC  methods  of 
arbitrary  order  can  be  constructed  that  retain  the  desired  stability  properties  necessary  for 
the  solution  of  stiff  differential  equations.  This  makes  it  an  appropriate  tool  for  the  solution 
of  parabolic  PDEs  in  time.  For  example,  recent  work  using  SDC  for  computational  fluid 
dynamics  include  [28]  [30]  EH- 

In  general,  implicit  methods  for  the  solution  of  a  system  of  linear  ODEs  yield  a  linear 
system  that  must  be  solved  on  each  iteration.  It  is  well  known  that  the  linear  system 
constructed  from  the  spatial  discretization  of  the  variable  coefficient  operator  in  (1.1)  is 
sparse.  In  M1,  this  linear  system  is  either  banded  or  block  banded  and  can  be  solved  in  0(n) 
where  n  is  the  number  of  spatial  discretization  nodes.  However,  in  M2,  the  straightforward 
spatial  discretization  leads  to  matrices  with  dimensionality  n2  x  n2  and  bandwidth  n.  While 
there  exist  fast  inversion  schemes  that  cost  0(n3),  alternating  direction  implicit  (ADI) 
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methods  replace  the  single  two-dimensional  implicit  step  with  two  sub-steps  where  only  one 
direction  is  treated  implicitly.  This  approach  results  in  schemes  with  computational  cost 
0(n2).  Likewise,  ADI  methods  in  M3  have  computational  cost  0(n3).  However,  popular 
ADI  methods  such  as  those  by  Peaceman-Rachford  or  Yanenko  are  low-order  (see,  e.g., 
[3911101111]).  In  this  dissertation,  we  combine  the  SDC  methods  for  the  solution  of  a  system 
of  ODEs  with  an  ADI  method  to  generate  computationally  efficient,  high-order  schemes  for 
the  solution  of  parabolic  PDEs  in  M2  and  M3. 


We  also  construct  high-order  representations  of  the  variable  coefficient  operator  in  ( 1.1 ). 


Traditionally,  pseudospectral  schemes  have  been  used  for  this  purpose;  however,  some  of 
the  eigenvalues  of  the  resulting  differentiation  operator  can  be  excessively  large.  An  im¬ 
provement  on  the  traditional  approach  is  to  subdivide  the  entire  spatial  domain,  construct 
bases  on  each  subdomain,  and  combine  the  obtained  discretization  with  the  implicit  SDC 
schemes.  The  resulting  schemes  are  high-order  and  have  CPU  time  requirements  that  are 


linear  in  the  number  of  spatial  discretization  nodes.  In  addition  to  solving  (1.1),  this  ap¬ 


proach  is  a  foundation  for  solving  nonlinear  equations  with  diffusion  components  as  well  as 
linear  and  nonlinear  problems  of  wave  propagation. 

This  dissertation  is  organized  as  follows.  Section  2  summarizes  various  standard  math¬ 
ematical  facts  to  be  used  in  subsequent  sections.  Section  3  contains  a  description  of  the 
SDC  methods  for  a  system  of  ODEs.  In  Section  4,  we  describe  low-order  solution  methods 
for  parabolic  PDEs  and  their  extensions  to  higher  dimensions.  In  Section  5,  we  discuss  the 
high-order  discretization  of  the  spatial  variable  in  order  to  construct  accurate  representa¬ 
tions  of  the  second-order  differentiation  matrix.  In  Section  6,  we  describe  the  SDC  methods 
for  parabolic  PDEs  and,  in  Section  7,  illustrate  the  behavior  of  these  methods  with  several 
numerical  examples. 

2  Preliminaries 

In  this  section,  we  introduce  notation  and  summarize  several  well  known  facts  which  we  use 
in  the  dissertation. 


2.1  Accuracy  and  Stability  of  ODE  Solvers 

We  assume  that  the  initial  value  problem  to  be  solved  is 

=  F(t,ip(t)),  t  G  [a,  b] 

ip(a)  =  Lpa 

where  <p(t)  :  M  — >  Cn  and  F  :  M  X  Cn  — >  Cn  is  sufficiently  smooth  (to  permit  high-order 
methods) . 


(2.1) 


Given  a  numerical  solution  to  (2.1),  <p(b),  the  numerical  method  is  said  to  be  of  order  k 
if,  for  any  sufficiently  smooth  F,  there  exists  a  constant  M  >  0,  such  that 

Mb)-v(b)\\<M(b-a)k+1. 

The  scalar  linear  differential  equation 

<//(t)  =  A t  >  0 

¥>(0)  =  1, 


(2.2) 

(2.3) 
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where  A  £  C,  has  exact  solution 


so  that  for  IZe(X)  <  0 


¥>(*)  =  eAt, 


lim  ip(t)  =  0. 

t->-oo 


(2.4) 


(2.5) 


Conceptually,  a  numerical  method  for  the  solution  of  (2.3)  with  IZe(X)  <  0  should  likewise 
yield  a  solution  that  converges  to  zero.  In  this  sense,  the  numerical  method  should  exhibit 
the  same  asymptotic  behavior  as  the  exact  solution  (see,  e.g. 


Any  numerical  method  for  the  solution  of  (2.3)  produces  a  sequence  of  approximations 
(pn  to  (p(nh )  for  n  =  0, 1, ... ,  which  satisfy  the  recurrence  relation 


(pn+i  =  g{hX)(pr, 


(2.6) 


Here,  h  either  denotes  the  step  size  for  multistep  schemes  or  the  interval  size  for,  e.g.,  Runge- 
Kutta  schemes.  Introducing  the  notation  z  =  hX,  the  stability  region  of  the  numerical 
method  is  the  set  of  all  z  E  C  such  that 


lim  ifn  =  0.  (2-7) 

n— >oo 

In  other  words,  the  stability  region  (together  with  its  boundary)  is  defined  to  be 

{zeC:\g(z)\  <1}  (2.8) 

where  g(z)  is  known  as  the  amplification  factor.  If  a  method  is  stable  for  all  z  in  the  left-half 
plane  (i.e. ,  IZe(z)  <  0),  then  it  is  said  to  be  A-stable.  A  method  is  said  to  be  A(a)-stable 
if  it  is  stable  for  all  z  such  that 


7T  —  a  <  arg (z)  <  ir  +  a.  (2-9) 

Thus,  A-stability  is  A(a)-stability  with  a  =  Finally,  if  a  method  is  A-stable,  it  is  said 
to  be  L-stable  if 

lim  g(z)  =  0.  (2.10) 

lZe(z)^—oo 

In  cases  when  g(z)  is  not  known  analytically,  the  stability  properties  of  the  method  can  be 
computed  numerically.  Specifically,  the  amplification  factor  g(z)  can  be  evaluated  as 

g(z)  =  <f(h).  (2.11) 


See  Section  I3i5l  for  further  details. 

For  a  given  e,  the  accuracy  region  is  defined  as  the  set  of  all  z  E  C  such  that 


eA  —  <p(h) 


<  £. 


Consider  now  the  system  of  linear  differential  equations 


ip'(t)  =  Aijj(t),  t  >  0 
V’(O)  =  T, 


(2.12) 


(2.13) 
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with 


(2.14) 


i>(t)  =  (V’i  , 


and 


1  =  {1,...,1}T.  (2.15) 

We  assume  that  A  is  a  diagonalizable  d  x  d-matrix  with  eigenvalues  Ai, . . . ,  Ad  such  that 
TZe(Xj)  <  0  for  all  j  =  1, . . .  ,d.  Then,  the  exact  solution  of  (2.13)  is 

< p{t )  =  eAt .  (2-16) 


Any  numerical  method  for  the  solution  of  (2.13)  produces  a  sequence  of  approximations  f>n 

(2.17) 


to  ip(nh )  for  n  =  0, 1, . . . ,  which  satisfy  the  recurrence  relation 

i>n+i  =  g{hA)%l)n. 


As  in  the  scalar  case,  g{z )  is  the  amplification  factor  and  the  definitions  for  A-stability, 
A(a)-stability,  and  L-stability  are  identical.  In  particular,  the  numerical  method  is  stable 
if  the  spectral  radius, 

p{g{hA))  <  1.  (2.18) 


2.2  Lagrange  Interpolation 

Given  nodes  to, ,  tm-i,  the  Lagrange  basis  polynomials 


satisfy  the  property 


Hence, 


m— 1 


ts  w = n 


i= 0 

*7 


t~tj 


j  =  0, . . . ,  m  -  1 


f'ji.ti)  —  Sij 


1  i  =  j 
0  i  /  j. 


m— 1 


3=0 


(2.19) 


(2.20) 


(2.21) 


is  the  Lagrange  interpolating  polynomial  for  f(t),  i.e. ,  ip(tj)  =  f(tj).  If  /  is  a  polynomial 
of  degree  m  —  1,  then  =  f(t). 


Observation  1.  For  evaluation  of  (2.21),  we  compute  £j(t)  as 


7/1  • 

m  = 


(2.22) 


where 


1 


and 


m—1 

0(*)=n  (*-**)•  (2-24) 

i=0 

iH 

Thus, 

m—1 

V’W  =  fi(t)  Y  ,  (2-25) 

3=0  3 

which  is  known  as  the  modified  form  of  Lagrange ’s  interpolation  formula.  The  computation 
of  each  Wj  costs  0(m),  so  {wj}  can  be  pre-computed  and  the  total  cost  is  0(m2)  operations. 
Once  { Wj }  have  been  computed,  the  cost  of  evaluating  fi(t)  at  any  point  t  is  0(m).  Since 
m  is  0(1),  the  cost  of  evaluating  the  function  f(t)  is  0(1). 


2.3  Legendre  Polynomials 

In  what  follows,  we  will  denote  by  Pj(x)  the  jth  Legendre  polynomial  on  the  interval  [—1,1], 
defined  by  the  three-term  recurrence 

pj+i{x)  =  2j  +  ixPj(x)  ~  ~jr[pi-i(x)  (2-26) 

with  the  initial  conditions  Pq(x)  =  1  and  Pi(x)  =  x  (see,  e.g.,  mm)-  Each  Legendre 
polynomial  satisfies  the  differential  equation 

.  r,d2Pj(x)  dPAx)  ,  .  „  ,  ,  , 

(1  ~  x)^y  ~  2x*r + j ' u  +  l)Pj{x)  =  °-  (2-27) 


The  family  of  polynomials  {-Pj}  is  orthogonal  on  [—1,1],  that  is, 

f1  2 

J  Pi(x)Pj(x)dx  =  — —5ij. 


We  denote  by  Pj  the  normalized  Legendre  polynomials,  that  is, 

Pj(x)  =  \l^4^pj(x)- 


(2.28) 


(2.29) 


The  following  Lemma  is  well  known.  Its  proof  can  be  found  in,  e.g.,  [3 j . 
Lemma  2.  Pj(x)  satisfies  the  following  two  relations: 


(2j  +  i)pfix)  =  p;+1(x)-p;_1(x) 


(2.30) 


and 

Pj(x)  =  (2 j  -  l)Pj-i(x)  +  (2 j  -  5 )Pj-3(x)  +  (2 j  -  9 )Pj-5(x)  +  . . .  (2.31) 


for  j  >  1. 
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2.4  Gaussian  Integration 

The  nodes  and  weights  of  Gaussian  quadratures  are  chosen  so  that  polynomials  of  degree 
less  than  or  equal  to  2 n  —  1,  {f o,  ■  •  • ,  are  integrated  exactly,  that  is, 


n— 1 


i=0 


(2.32) 


for  all  j  =  0,  ...,2n  —  1.  Provided  that  the  function  /  is  well  approximated  by  linear 
combinations  of  such  polynomials,  the  Gaussian  quadrature  formulae  provide  good  approx¬ 
imations  to  integrals  of  the  form 

f(t)oj(t)dx,  (2.33) 

where  oj(x)  is  an  integrable  non-negative  function  mm- 

The  following  are  common  choices  of  quadrature  nodes  to, . . .  ,tn-i  for  the  integration 
of  functions  defined  on  [a,b]  =  [—1, 1].  In  all  cases,  ui(t)  =  1,  except  for  Chebyshev  nodes 
for  which  u(t)  =  (1  —  t2)“T  See,  e.g.,  [H  IS  ;8l  [39]  for  further  details. 

•  Gauss-Legendre:  ti  is  the  Ith  zero  of  Pn(t ),  the  nth-degree  Legendre  polynomial. 

•  Chebyshev:  ti  =  cos(^p7r),  which  is  the  ith  zero  of  Tn(t ),  the  nth- degree  Chebyshev 
polynomial. 

•  Gauss-Lobatto:  to  =  —1,  in-i  =  1 ,  and  ti  is  the  ith  zero  of  P^_i(t)  for  i  =  2, . . . ,  m  —  2. 

•  Left-hand  Gauss- Radau:  to  =  —1  and  ti  is  the  ith  zero  of 

Pn_i(t)  +  Pn(t) 

1  +  t  v  ' 

The  right-hand  Gauss-Radau  quadratures  are  the  negative  of  the  left-hand  Gauss- 
Radau  quadratures. 


In  this  dissertation,  we  will  also  use  what  we  call  a  right-hand  variant  of  Gauss-Legendre 
nodes, 


*■  -  2±tti  - 1  (2'35) 

where  yi  is  the  ith  zero  of  Pn{t).  We  refer  to  it  as  a  right-hand  variant  of  Gauss-Legendre 
because  xn-\  =  1. 

Observation  3.  In  some  cases,  analytical  expression  for  the  weights  Wi  exist.  However,  in 
general,  once  the  nodes  to,  ■  ■  ■  ,tn-i  are  computed,  the  weights  wq,  . . . ,  wn-i  can  be  deter¬ 
mined  by  solving  the  linear  system  of  equations 

n— 1  -l 

J2wiMti)=  /  <fj(t)uj(t)dt 
i= o  J~l 

for  i  =  0, . . . ,  n  —  1. 
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3  Spectral  Deferred  Corrections 

3.1  The  Picard  Integral  Equation 

Integrating  the  initial  value  problem  in  0  with  respect  to  t,  yields  the  Picard  integral 
equation 


(p{t)  =tpa+  F(r,^(r))dr. 


(3.1) 


Suppose  we  have  an  approximate  solution  (p{t )  to  (2.1).  A  measure  of  the  quality  of  the 
approximation  is  given  by  the  residual  function 


e(t)  =  pa  +  F(t,  <p(r))dT  -  (pit). 


We  define  the  error  5(t)  by 


S(t)  =  Pit )  -  pit). 


Substituting  (3.3)  back  into  (3.1),  we  obtain 


pit)  +  5it)  =  pa+  F(r,  (^(r)  +  d(r))dr. 


Substituting  (3.2)  into  (3.4)  yields 

S(t)  =  f  [Fir,  p(t)  +  <5(t))  -  F(t,  <£(t))]  dr  +  e(t). 

J  a 

By  defining  the  function  G  :  M  x  Cn  — >  Cn  to  be 

G{t,  <5)  =  F{t,  pit)  +  5(f))  -  Fit,  pit)), 


we  can  rewrite  (3.5)  as  a  Picard-type  integral  equation 


(5(f)—  /  G(r,  (5(r))dr  =  e(f). 


(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 


This  is  the  key  equation  for  Spectral  Deferred  Correction  (SDC)  methods  as  in,  e.g.  [TU1I2?]. 

3.2  Implicit  Euler  for  the  Picard  Equation 

Suppose  that  we  have  nodes  a  <  to  <  ti  <■■■  <  tm-\  <  b.  Then,  the  implicit  Euler  scheme 


(3.8) 


for  the  solution  of  (2.1),  or  equivalently  (|3.l[),  is  given  by 

Pi+l  —  Pi  T  hi  •  Fiti+\,Pi+\),  hi  —  ti+l  ti 


for  i  =  0, . . . ,  m  —  2. 

Similarly,  the  implicit  Euler  scheme  for  the  solution  of  (3.7)  is  given  by 

^i+i  =  Si  +  hi  ■  G(tj_|_ i,  Pi+i)  +  (^(tj-i-i)  —  e(tj)).  (3-9) 

Observation  4.  Given  an  approximate  solution  (pit )  to  \2.1\),  a  key  component  of  the 


scheme  in  (3.9)  is  the  evaluation  of  the  integral  in  (3.2).  This  requires  stable,  high-order 
methods  for  interpolation  and  integration. 


3.3  Spectral  Integration 

We  would  like  to  evaluate  the  integral 


(3.10) 


for  i  =  0, . . . ,  m  —  1.  A  standard  approach  is  to  represent  the  function  /  by  a  truncated 
series  expansion  and  compute  the  integral  by  analytically  integrating  the  series.  The  result 
is  a  linear  mapping  between  the  vector  of  m  function  values  and  the  vect°r 

ft  1  m—X 

l  fa  /(r)^T  f  •  This  approach  is  numerically  stable  provided  that  the  quadrature  nodes 

are  selected  carefully,  e.g.,  Gaussian-Legendre  or  Chebyshev  nodes  (see  Section  2.4).  In  this 
case,  the  matrix  representation  of  the  linear  mapping  is  referred  to  as  the  spectral  integration 
matrix.  Its  singular  values  are  between  0(1), . . . ,  0(l/m2).  See,  e.g.,  ([UlllHl  39])  for  a 
more  detailed  discussion  on  spectral  integration  matrices  and  their  numerical  properties. 

The  classical  approach  to  the  construction  of  the  spectral  integration  matrix  is  to  use 
Lagrange  interpolating  polynomials  defined  in  (2.21).  Specifically,  given  nodes  to, . . .  ,fm_i 
and  a  function  /(f),  the  integral  of  the  Lagrange  interpolating  polynomial, 


m—  1 

V>(t)  =  (3-n) 

3=0 


is 


pti  m—  1 

/  if{r)dT  =  f(tj)  tj(r)dT. 

J  a  j_  o  J  cl 

(3.12) 

If  we  define 

II 

's) 

(3.13) 

then 

rti 

/  V’C T)dr  =  Sf 

J  a 

(3.14) 

where  /  is  a  vector  with  elements  /,;  =  /(tj).  The  matrix  S  is  referred  to  as  the  integration 


matrix.  If  the  nodes  to, . . . ,  fm_i  are  Gauss-Legendre  or  Chebyshev  nodes,  then  S  is  referred 
to  as  the  spectral  integration  matrix. 


Observation  5.  An  alternative  formulation  useful  for  the  computation  of  the  spectral  in¬ 
tegration  matrix  in  (3.13)  expresses  the  function  f(t )  through  Legendre  polynomials  instead 

so  that 


of  Lagrange  interpolating  polynomials. 


Specifically,  given  {aj}JL f 


m—  1 


f(t)  =  a-jPjft), 
3=0 


we  woidd  like  to  compute  {/3j}jL0  so  that 


(3.15) 


(3.16) 
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Denoting 


we  have 
where 


f  =  if  (to),  •  •  • 

«  =  {«()>  ■  •  •  ,  «m—  l}T  , 

7  =  Va 


Vij  =  Pjiti). 

To  compute  the  indefinite  integral,  we  have 

r-t  m~1  rt 


/t  rt 

f  (t) dr  =  /  pj(T)dr. 

1  j=o  -7-1 


Using  (2.30),  we  obtain 

rn-l  rt 


m— 1 


Y  ai  /  Pi(T)dT  =  a°  (Po(i)  +  Pi(t))  +  Y  T~rj  (Pi+1(t)  “  Pi- 1(*))  • 

j=o  J-1  j=i  J  + 


Hence,  we  arrive  at 


Pi  =  ' 


ao  —  |ai  j  =  0 


Otj- 1  +  1 

2j— 1  2j+3 

Otm  —  2 

2m— 3 
C^m  —  1 


0  <  j  <  m  —  1 
j  =  m  —  1 
j  =  m. 


_  2m— 1 

The  denote  by  A,  the  {m  +  1)  X  m-matrix  such  that 

Aa  =  (3, 

where 

P  =  {Poi  ■  ■  ■ ,  Pm}T  ■ 

Finally,  we  denote  by  V  the  m  x  (m  +  1 ) -matrix  such  that 

Vi;j  =  Pj{ti). 


(3.17) 

(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 

(3.23) 

(3.24) 

(3.25) 


Then,  it  is  straightforward  to  see  that  the  spectral  integration  matrix,  S,  is  given  by 

S  =  VAV-1.  (3.26) 

Given  any  approximate  solution  (p  (f),  the  residual  function  is 

e(t)  =  ipa+[  F(t,  <p(r))dT  -  (pit).  (3.27) 

J  a 

Denoting 

V  =  M^o),  •  •  -,<p(tm- l)}T 

Tp-a  =  {pa,...,ipa}T  (3.28) 

F(<P)  =  {p{t,  v(to)),  ■  ■  •  ,  F(tm- 1,  ^(tm-l))}T  , 

we  approximate  the  residual  function  by  the  vector 

aifi)  =  tpa  +  SF(ip)  -  Tp  (3.29) 

where  S  is  the  spectral  integration  matrix. 
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3.4  Spectral  Deferred  Corrections  Algorithm 

Spectral  Deferred  Correction  methods  for  the  solution  of  (|2.1[)  proceed  as  follows.  First,  an 
approximate  solution  ip°{t)  with  error  6(t )  =  ip(t)  —  ip°(t)  is  computed  at  quadrature  nodes 
a  <  to  <  t\  <■■■  <  tm-i  <  b  via  some  numerical  method  -  generally,  low  order.  Then, 
an  approximation  for  the  error  8(t)  is  computed  at  the  quadrature  nodes  and  the  updated 
solution  becomes  ipl{t)  =  ip°{t)  +  5(f).  This  procedure  is  repeated  L  times  to  provide  a 
more  accurate  approximation  ipL(t). 

More  explicitly,  given  some  choice  of  quadrature  nodes,  Algorithm [l]  is  the  implicit  SDC 
algorithm  from  [10]. 


Algorithm  1  Implicit  Spectral  Deferred  Corrections 


Compute  an  initial  approximate  solution,  <p°(tj ),  for  j  =  0, . . . ,  m  —  1  via  the  implicit  Euler 


method  (3.8). 


For  l  =  1 , ...  ,L 


1.  Compute  the  residual  function  e(tj)  in  (3.2)  via  spectral  integration  (see  Section  3.3), 


2.  Compute  5(tj)  via  the  implicit  Euler  scheme  (3.9), 

3.  Update  the  approximate  solution  if J(tj )  =  c pl~l(tj)  +  5(tj). 


Set  the  resulting  solution  as  <p(tj)  =  ipL(tj). 


Given  the  value  of  the  solution  ip(t)  at  to, . . .  ,tm_ i,  we  can  evaluate  the  solution  every¬ 
where  on  the  interval  [a,  b]  via  Lagrange  interpolation  (see  Section  2.2). 


Observation  6.  Many  choices  of  to, . . . ,  tTO-i  do  not  include  the  right  endpoint.  If  (2.1) 
is  not  stiff,  then  ip(b)  can  be  computed  via  Gaussian  integration  (see  Section  2-4)  as 


m—  1 


ip(b)  =ipa+  F(t ,  ip{r))dT  =  ipa  +  ^2  WiF(ti ,  ip(ti)). 


(3.30) 


i= 0 


For  stiff  systems,  this  approach  is  numerically  unstable.  Instead,  the  solution  can  be  com¬ 
puted  at  t  =  b  via  Lagrange  interpolation.  Nevertheless,  it  is  for  this  reason  that  quadrature 
schemes  that  include  the  right  end  point  (i.e.,  tm- 1  =  b)  are  preferred  (e.g.  right-hand 
Gauss-Radau  or  right-hand  Gauss-Legendre) . 

See,  e.g.,  mm  for  the  proof  of  the  following  theorem. 

Theorem  7.  For  any  sufficiently  smooth  function  F,  the  approximate  solution  computed 
via  SDC  converges  to  the  exact  solution  with  order  min(m,  L  +  1). 

Observation  8.  Replacing  the  implicit  Euler  method  with  the  higher-order  trapezoidal 
method  reduces  the  required  number  of  corrections,  but  results  in  a  scheme  that  is  not  L- 
stable.  Moreover,  such  an  approach  requires  careful  implementation.  For  example,  using 
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Gauss- Legendre  nodes  and  extrapolating  the  solution  to  t  =  b  via  Lagrange  interpolation  can 
result  in  a  scheme  whose  stability  properties  are  like  that  of  an  explicit  method.  We  avoid 
such  concerns  by  using  the  implicit  Euler  method  to  drive  the  SDC  algorithm. 


3.5  Stability  and  Accuracy  Regions  of  Spectral  Deferred  Corrections 


In  [TU],  Spectral  Deferred  Corrections  using  Gauss-Legendre  quadratures  nodes  are  intro¬ 
duced.  However,  the  choice  of  the  nodes  strongly  influences  the  stability  and  accuracy 
regions  of  the  resulting  schemes.  In  this  section,  we  compare  the  stability  and  accuracy 
properties  of  SDC  methods  with  four  classes  of  nodes:  Gauss-Legendre,  Chebyshev,  Gauss- 


Lobatto,  and  what  we  call  a  right-hand  variant  of  Gauss-Legendre  (see  Section  2.4). 

we  introduce  the  concepts  of  stability  and  accuracy  regions  of  numerical 


In  Section  2.1 


methods  for  the  solution  of  ODEs.  For  clarity,  we  restate  several  facts  here.  To  compare  the 
stability  and  accuracy  properties  of  the  SDC  methods,  we  apply  them  to  the  scalar  linear 
differential  equation 

<p'(t)  =  A ip(t),  t>  0 

¥>(0)  =  1, 


where  A  E  C,  has  exact  solution 


¥>(*)  =  e 


A  t 


(3.32) 


Traditionally,  for  a  fixed  time  step  h,  the  stability  of  the  numerical  method  is  formulated 


in  terms  of  the  so-called  amplification  factor  g(hX)  in  (2.6).  The  stability  region  (together 
with  its  boundary)  is  defined  to  be 


{z  G  C  :  \g{z)\  <  1}.  (3.33) 

where  z  =  hX.  If  a  method  is  stable  for  all  z  such  that  IZe(z)  <  0,  then  it  is  said  to  be 
A-stable.  If  the  method  is  A-stable,  it  is  said  to  be  L-stable  if 


lim  g(z)  =  0.  (3.34) 

lZe(z)^—oo 

Finally,  for  a  given  e,  the  accuracy  region  is  defined  as  the  set  of  all  z  E  C  such  that 


e 


A 


<p(h) 


<  £, 


(3.35) 


where  <p(h)  is  the  approximate  solution  at  t  =  h  computed  via  some  numerical  method. 
Since  the  stability  function  of  SDC  methods  is  not  known  analytically,  the  stability  region 
of  the  methods  are  computed  numerically  via  Algorithm  [2j  A  simple  modification  of  the 
algorithm  is  used  to  compute  the  accuracy  region. 

If  we  fix  h  =  1,  then  a  simple  calculation  shows  that  for  a  given  A,  (p{\)  =  g( A)  =  g(z). 
Therefore,  to  compute  the  boundary  of  the  stability  region,  we  find  A  such  that  ip(  1)  =  1. 
To  start,  we  use  the  fact  that  for  z  =  0,  g(z)  =  1  for  SDC  methods.  The  algorithm  then 
proceeds  to  take  small  steps  and,  after  each  step,  searches  in  the  direction  perpendicular  to 
that  step  to  find  A  such  that  cp(l)  =  1. 
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Algorithm  2  Stability  Region  Boundary  Search 

Initialization: 

1.  Starting  from  Ao  =  0,  take  a  small  step  along  the  positive  imaginary  axis  and  denote 
this  point  Ai, 

2.  For  A  =  Ai,  determine  if  (p(l)  >  1  or  ip(  1)  <  1, 

3.  If  ip(  1)  >  1,  search  parallel  to  the  real  axis  in  the  negative  direction  to  find  A  such 
that  <f(l)  =  1.  If  v3(l)  <  1,  search  parallel  to  the  real  axis  in  the  positive  direction  to 
find  A  such  that  (p{  1)  =  1.  Denote  this  by  Ai. 

For  j  =  1,2, . . . 

1.  Starting  from  A  j,  take  a  small  step  in  the  same  direction  as  the  vector  connecting 
Aj_i  to  A  j  and  denote  this  point  A  j.  Denote  the  vector  by  Vj  =  { Xj,yj}T , 

2.  For  A  =  Xj,  determine  if  (p{  1)  >  1  or  </?(  1)  <  1, 

3.  If  ip(  1)  >  1,  search  in  the  direction  {— yj,  Xj}T  to  hnd  A  such  that  (p(l)  =  1.  If 
(p{  1)  <  1,  search  in  the  direction  {yj,  —Xj}T  to  find  A  such  that  <p(  1)  =  1.  Denote  this 
by  A  j. 

Continue  the  iteration  until  A  returns  to  the  origin. 


In  this  section,  we  compare  SDC  methods  of  order  4,  6,  8,  10,  and  12.  In  each  case, 
the  number  of  nodes  m  equals  to  the  order  of  the  SDC  method  and  the  number  of  deferred 
corrections  L  equals  to  m  —  1.  Similar  analysis  for  other  choices  of  nodes  can  be  found  in 
[29].  For  a  given  order,  we  compare  the  stability  region,  the  amplification  factor  g(z)  as 
7 Ze(z)  — >  — oo,  and  the  accuracy  region  for  different  discretization  nodes.  Unless  otherwise 
stated,  in  the  plots  of  the  stability  region  in  Figures  3.1  3.4  3.7  |3.10  and  3.13  the  lines 
indicate  where  g(z)  =  1  with  g(z)  <  1  outside  of  the  lines.  In  the  plots  of  the  amplification 
factor  g(z)  for  real  z  in  Figures  3.2,  3.5,  3.8,  3.11,  and  3.14  x  =  —  log10(— z).  Thus,  the 
x-axis  in  these  plots,  which  has  range  x  G  [—8, 5],  corresponds  to  z  G  [ —  10s ,  — 10-5].  Finally, 
in  the  plots  of  the  accuracy  region  in  Figures  3.  3j  |3.6|  |3.9]  |3.12]  and  |3.15]  the  lines  indicate 
where  e  =  1CD5. 


In  Figures  3.1,  3.2 


and  |3.3|  we  compare  the  stability  and  accuracy  properties  of  SDC 
methods  of  order  4  for  different  sets  of  discretization  nodes.  Figure  |3.1|  is  a  plot  of  the 
stability  regions,  Figure  3.2  is  a  plot  of  the  amplification  factors  g(z)  as  7 Ze(z)  —  oo,  and 
Figure [3~3] is  a  plot  of  the  accuracy  regions.  In  this  case,  the  fourth  order  SDC  method  with 
Gauss-Legendre  nodes  is  A-stable,  while  the  rest  are  A(a)-stable  with  a  slightly  less  than 
|.  For  all  choices  of  discretization  nodes,  g(z)  <  1  as  7 Ze(z)  — >•  — oo.  However,  clearly  the 
SDC  method  with  Gauss-Lobatto  nodes  has  a  much  larger  amplification  factor  in  the  limit 
than  the  other  choices  of  nodes. 


13 


-0.002  -0.0015  -0.001  -0.0005  0  0.0005  0.001  0.0015  0.002 


Figure  3.1:  Comparison  of  the  stability  regions  for  fourth  order  SDC  methods  with  different 
sets  of  discretization  nodes.  The  lines  indicate  where  the  amplification  factor  g(z)  =  1  with 
g(z)  <  1  outside  of  the  lines.  The  bottom  figure  is  a  plot  of  the  stability  regions  zoomed  in 
on  the  origin. 
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Figure  3.2:  Comparison  of  the  amplification  factor  g(z)  as  TZe(z)  — >  — oo  for  fourth  order 
SDC  methods  with  different  sets  of  discretization  nodes.  In  the  x-axis  of  the  plot,  x  = 

-logio(-z). 


Figure  3.3:  Comparison  of  the  accuracy  regions  with  e  =  10  5  for  fourth  order  SDC  methods 
with  different  sets  of  discretization  nodes. 


In  Figures  3.4,  3.5 


and  |3.6|  we  compare  the  stability  and  accuracy  properties  of  SDC 
methods  of  order  6  for  different  sets  of  discretization  nodes.  Figure  |3.4|  is  a  plot  of  the 
stability  regions,  Figure  3.5  is  a  plot  of  the  amplification  factors  g(z)  as  IZe(z)  — >  — oo, 
and  Figure  3.6  is  a  plot  of  the  accuracy  regions.  In  this  case,  none  of  the  sixth  order  SDC 
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methods  are  H-stable,  but  they  all  are  H(a)-stable  with  a.  slightly  less  than  For  all 
choices  of  discretization  nodes,  g(z)  <  1  as  IZe(z)  — >  —  oo.  However,  as  in  the  fourth  order 
case,  clearly  the  SDC  method  with  Gauss-Lobatto  nodes  has  a  much  larger  amplification 
factor  in  the  limit  than  the  other  choices  of  nodes. 


Figure  3.4:  Comparison  of  the  stability  regions  for  sixth  order  SDC  methods  with  different 
sets  of  discretization  nodes.  The  lines  indicate  where  the  amplification  factor  g(z )  =  1  with 
g(z)  <  1  outside  of  the  lines.  The  bottom  figure  is  a  plot  of  the  stability  regions  zoomed  in 
on  the  origin. 
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Figure  3.5:  Comparison  of  the  amplification  factor  g(z)  as  7 Ze(z)  —>■  — oo  for  sixth  order 
SDC  methods  with  different  sets  of  discretization  nodes.  In  the  x-axis  of  the  plot,  x  = 
l°gio(  z) ■ 


Figure  3.6:  Comparison  of  the  accuracy  regions  with  e  =  10  5  for  sixth  order  SDC  methods 
with  different  sets  of  discretization  nodes. 


In  Figures  3.7,  3. 


and  |3.9|  we  compare  the  stability  and  accuracy  properties  of  SDC 
methods  of  order  8  for  different  sets  of  discretization  nodes.  Figure  |3.7|  is  a  plot  of  the 
stability  regions,  Figure  3.8  is  a  plot  of  the  amplification  factors  g(z)  as  TZe(z)  — >  — oo,  and 
Figure  |3.9|  is  a  plot  of  the  accuracy  regions.  In  this  case,  none  of  the  eighth  order  SDC 
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methods  are  H-stable,  but  they  all  are  H(a)-stable  with  a  slightly  less  than  For  all 
choices  of  discretization  nodes,  g(z )  <  1  as  TZe(z)  — >  —  oo.  However,  as  in  the  fourth  and 
sixth  order  cases,  clearly  the  SDC  method  with  Gauss-Lobatto  nodes  has  a  much  larger 
amplification  factor  in  the  limit  than  the  other  choices  of  nodes. 


- 

Legendre  - 

Chebyshev  - 

Right-hand  Legendre  . 

Lobatto 

- 

- 

- 

- 

-0.002  -0.0015  -0.001  -0.0005  0  0.0005  0.001  0.0015  0.002 


Figure  3.7:  Comparison  of  the  stability  regions  for  eighth  order  SDC  methods  with  different 
sets  of  discretization  nodes.  The  lines  indicate  where  the  amplification  factor  g(z)  =  1  with 
g(z)  <  1  outside  of  the  lines.  The  bottom  figure  is  a  plot  of  the  stability  regions  zoomed  in 
on  the  origin. 
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Figure  3.8:  Comparison  of  the  amplification  factor  g{z )  as  IZe(z)  — >  — oo  for  eighth  order 
SDC  methods  with  different  sets  of  discretization  nodes.  In  the  x-axis  of  the  plot,  x  = 
-logio(-z). 


Figure  3.9:  Comparison  of  the  accuracy  regions  with  e  =  10  5  for  eighth  order  SDC  methods 
with  different  sets  of  discretization  nodes. 


In  Figures [3. 10t  |3-lll  and|3.12|we  compare  the  stability  and  accuracy  properties  of  SDC 
methods  of  order  10  for  different  sets  of  discretization  nodes.  Figure  3.10|  is  a  plot  of  the 
stability  regions,  Figure  3.11  is  a  plot  of  the  amplification  factors  g(z)  as  7 Ze(z)  — >  —  oo,  and 
Figure  |3.12  is  a  plot  of  the  accuracy  regions.  In  the  plot  of  the  stability  regions  in  Figure 
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3.10,  the  lines  indicate  where  g(z)  =  1  with  g(z)  <  1  outside  of  the  lines  except  for  the  SDC 
method  with  Gauss-Lobatto  nodes  where  the  line  indicated  where  g(z )  =  1  with  g(z)  <  1 
inside  of  the  line.  In  this  case,  none  of  the  tenth  order  SDC  methods  are  ^4-stable,  but, 
with  the  exception  of  the  SDC  method  with  Gauss-Lobatto  nodes,  they  all  are  yf(a)-stable 
with  a  slightly  less  than  Despite  the  implicit  construction  for  the  SDC  scheme  with 
Gauss-Lobatto  nodes,  the  stability  region  is  like  that  of  an  explicit  method.  For  all  choices 
of  discretization  nodes  except  Gauss-Lobatto,  g(z )  <  1  as  IZe(z)  -*  — oo. 
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Figure  3.10:  Comparison  of  the  stability  regions  for  tenth  order  SDC  methods  with  different 
sets  of  discretization  nodes.  The  lines  indicate  where  the  amplification  factor  g(z)  =  1  with 
g(z)  <  1  outside  of  the  lines  except  for  Gauss-Lobatto  nodes  for  which  g(z)  <  1  inside  of 
the  line.  The  bottom  figure  is  a  plot  of  the  stability  regions  zoomed  in  on  the  origin. 
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Figure  3.11:  Comparison  of  the  amplification  factor  g(z)  as  1Ze(z )  — >  — oo  for  tenth  order 
SDC  methods  with  different  sets  of  discretization  nodes.  In  the  x-axis  of  the  plot,  x  = 
l°gio(  z). 


Figure  3.12:  Comparison  of  the  accuracy  regions  with  e  =  10  5  for  tenth  order  SDC  methods 
with  different  sets  of  discretization  nodes. 


In  Figures  [3. 13|  |3. 14[  and |3.15|  we  compare  the  stability  and  accuracy  properties  of  SDC 
methods  of  order  12  for  different  sets  of  discretization  nodes.  Figure  3.13|  is  a  plot  of  the 
stability  regions,  Figure  3.14  is  a  plot  of  the  amplification  factors  g(z )  as  IZe(z)  — >  —  oo,  and 
Figure  |3.15  is  a  plot  of  the  accuracy  regions.  In  the  plot  of  the  stability  regions  in  Figure 
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3.13,  the  lines  indicate  where  g(z)  =  1  with  g(z)  <  1  outside  of  the  lines  except  for  the  SDC 
method  with  Gauss-Lobatto  nodes  where  the  line  indicated  where  g(z )  =  1  with  g(z)  <  1 
inside  of  the  line.  In  this  case,  none  of  the  twelfth  order  SDC  methods  are  A-stable,  but, 
with  the  exception  of  the  SDC  method  with  Gauss-Lobatto  nodes,  they  all  are  yf(a)-stable 
with  a  slightly  less  than  Despite  the  implicit  construction  for  the  SDC  scheme  with 
Gauss-Lobatto  nodes,  the  stability  region  is  like  that  of  an  explicit  method.  For  all  choices 
of  discretization  nodes  except  Gauss-Lobatto,  g(z )  <  1  as  IZe(z)  -*  — oo. 


- 

Legendre  - 

Chebyshev  - 

Right-hand  Legendre  . 

- 

■ 

- 

- 

- 

-0.002  -0.0015  -0.001  -0.0005  0  0.0005  0.001  0.0015  0.002 


Figure  3.13:  Comparison  of  the  stability  regions  for  twelfth  order  SDC  methods  with  differ¬ 
ent  sets  of  discretization  nodes.  The  lines  indicate  where  the  amplification  factor  g(z)  =  1 
with  g{z)  <  1  outside  of  the  lines  except  for  Gauss-Lobatto  nodes  for  which  g(z)  <  1  inside 
of  the  line.  The  bottom  figure  is  a  plot  of  the  stability  regions  zoomed  in  on  the  origin. 
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Figure  3.14:  Comparison  of  the  amplification  factor  g{z )  as  IZe(z)  — >  — oo  for  twelfth 
order  SDC  methods  with  different  sets  of  discretization  nodes.  In  the  x-axis  of  the  plot, 
*  =  k>6io(  ~)- 


Figure  3.15:  Comparison  of  the  accuracy  regions  with  e  =  10  5  for  twelfth  order  SDC 
methods  with  different  sets  of  discretization  nodes. 


It  is  clear  from  the  stability  regions  in  Figures  |3.10|  and  3.13  that  the  choice  of  Gauss- 
Lobatto  nodes  for  high-order  schemes  is  not  a  practical  option.  Between  the  other  choices 
of  nodes,  the  differences  in  the  stability  and  accuracy  regions  are  fairly  small.  However, 
the  SDC  scheme  with  the  right-hand  variant  of  Gauss-Legendre  nodes  consistently  had 
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(slightly)  larger  accuracy  regions  with  essentially  equivalent  stability  properties.  For  this 
reason,  in  our  implementation  of  SDC  methods  for  evolving  parabolic  differential  equations 
in  time,  we  use  the  right-hand  variant  of  Gauss-Legendre  nodes. 


4  Parabolic  PDE  Solvers 


In  order  to  develop  an  SDC  algorithm  for  parabolic  PDEs,  we  need  a  low-order  numerical 
solver  to  drive  the  PDE  version  of  the  SDC  method.  This  scheme  serves  the  same  role 


as  that  the  implicit  Euler  method  in  (3.8)  and  (3.9)  does  for  the  implicit  SDC  algorithm 
for  the  solution  of  ODEs.  In  many  situations,  the  scheme  should  retain  the  L-stability 
property  provided  by  the  implicit  Euler  method  for  ODEs.  In  this  section,  we  discuss  the 
implicit  Euler  method  for  parabolic  PDEs  in  one  spatial  variable  and  its  extension  to  higher 
dimensions,  known  as  the  Alternating  Direction  Implicit  (ADI)  method.  Using  ADI  helps 
control  the  computational  cost  of  the  numerical  solver.  As  a  general  reference,  we  point  to, 
e.g.,  [321 1401  (33j  for  a  more  detailed  discussion  of  numerical  techniques  for  the  solution  of 
PDEs. 


4.1  Implicit  Euler  Methods  for  Parabolic  PDEs 

Consider  the  parabolic  PDE  in  one  dimension 


d(j) 

dt 


(4.1) 


where  a  is  the  diffusion  coefficient.  Discretizing  (4.1 )  using  spatial  nodes  xo, . . . ,  xn-\  yields 
a  system  of  ODEs 

^  =  DX  ( ADX )  (4.2) 

where 

<P  =  =  {4>(x0,  t),...,  4>(xn-1,t)}T  .  (4.3) 

Matrices  Dx  and  Dx  are  discrete  approximations  of  the  derivative  operator  and  A  is  the 
operator  of  point- wise  multiplication  of  the  function  a(x)  defined  via  the  formula 

A((p)(Xj)  =  a(xj)<p(xj).  (4.4) 

Observation  9.  Our  construction  of  Dx  and  Dx  results  in  a  matrix  DXADX  that  is  sym¬ 
metric.  While  the  matrices  Dx  and  Dx  both  approximate  the  derivative  operator,  in  our 
construction  they  are  not  necessarily  the  same.  In  particular,  the  matrix  Dx  enforces  the 
relevant  boundary  conditions,  while  the  matrix  Dx  does  not.  For  details,  see  Sections 
and\5.3i 

Given  nodes  to  <•••  <  tm- 1,  the  implicit  Euler  method  for  the  solution  of  (4.2)  is 
’(<+1)  =  +  hibxADx^i+1\  hi  =  ti+ 1  -  u 


5.2 


(4.5) 


or 


<p{i+1)  =  (I-  hiDxADxy\ (i).  (4.6) 

The  following  Lemma  is  an  immediate  consequence  of  the  definition  of  stability  in  ( 2. 18|) 
and  L-stability  in  (2.10). 
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Lemma  10.  The  implicit  Euler  method  is  stable  provided 


p  ((/  -  hiDxADx)-^  <  1  (4.7) 

where  p(-)  denotes  the  spectral  radius  of  a  matrix.  Hence,  it  is  stable  if  DXADX  is  negative 
semi- definite.  Finally,  the  implicit  Euler  method  is  also  L-stable. 

Observation  11.  For  most  practical  implementations  of  the  implicit  Euler  method  for 
parabolic  PDEs,  the  matrix  DXADX  is  banded  or  block  banded.  As  as  result,  it  is  possible  to 
efficiently  solve  linear  systems  with  the  matrix  I —hiDxADx.  In  this  case,  the  computational 
cost  of  solving  such  linear  systems  is  0(n).  For  parabolic  PDEs  in  M2,  the  discretized  spatial 
grid  is  nxn  and,  therefore,  the  matrix  representation  of  the  Laplacian  is  n2  xn2  -  see  Figure 


However,  we  instead  use  an  alternating  direction  implicit  scheme  (see,  e.g.,  l3fk\4U(),  which 
is  much  simpler  to  implement  and  more  computationally  efficient.  Using  an  alternating 
direction  implicit  scheme  is  even  more  advantageous  in  M3. 


4-1  In  this  case,  there  exist  fast  inversion  schemes  that  cost  0(n3)  (see,  e.g.,  |B  122,  HlFj- 


Figure  4.1:  Matrix  representation  of  the  Laplacian  in  M2. 


4.2  Alternating  Direction  Implicit  Methods  for 
Parabolic  PDEs 

Consider  the  parabolic  PDE  in  two  dimensions 


^  =  V-(a(s,j/)V$. 


(4.8) 


Discretizing  (|4.8|)  using  spatial  nodes  xo,  •  •  • ,  xn-i  and  yo, ... ,  yn-i  yields  a  system  of  ODEs 

(4.9) 
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—  (  Dx  ( ADX )  +  Dy  ( ADy ) )  ip, 


where 


V  =  <p(t)  =  <!>{xjx,yjy,t) 


(4.10) 


for  all  jx  =  0, ...  ,7i  —  1  and  jy  =  0, . . . , n  —  1.  The  n2  x  n2-matrices  DXADX  and  DyADy 
are  discrete  approximations  of  the  second-order  differentiation  operator  with  variable  coef¬ 
ficients.  For  details,  see  Sections  4.1|  5.2  and  5.3  In  (4.9),  A  is  the  operator  of  point-wise 
multiplication  of  the  function  a(x,  y )  defined  via  the  formula 


A(v)  (xja ,  yjy )  =  a(xjx ,  yjy)(p(xjx ,  yjy 


(4.11) 


Conceptually,  ADI  schemes  in  M2  replace  the  single  two-dimensional  implicit  step  with 
two  sub-steps  where  only  one  direction  is  treated  implicitly.  This  approach  replaces  a  single 
n2  x  n2  system  with  n  systems  of  size  nxnon  each  sub-step.  Since  each  system  can  be  solved 
with  computational  cost  0(n)  (see  Observation  11),  the  total  computational  cost  of  ADI 
schemes  in  M2  are  0(n2),  with  a  small  constant.  Likewise,  in  M3,  the  total  computational 
cost  is  0(n3). 

Given  nodes  t o,  •  •  • ,  tm-i ,  a  popular  ADI  method  is  the  so-called  Peaceman-Rachford 
ADI  method  (see,  e.g., 


(4.12) 


=  ptt  +  *±bxADx<p«+h)  +  ^DyADyifW 


<p(i+ 1)  =  +  -jDxADx^i+1^  +  ^DyADy^i+1\ 


We  rewrite  (4.12)  in  the  form 


where 


=L^\ 

(4.13) 

^DxADx)~1(I+^DxADx)  o 

A  A 

^DyADj-^+^DyADy). 

(4.14) 

The  following  Lemma  is  an  immediate  consequence  of  the  definition  of  stability  in  (2.18) 
and  L-stability  in  (|2. 10|) .  See,  e.g.,  [39l  1401 144]  for  further  details. 


Lemma  12.  The  Peaceman-Rachford  ADI  method  is  stable  provided 


p(L)  <  1.  (4.15) 

Hence,  it  is  stable  if  is  stable  if  both  DXADX  and  DyADy  are  negative  semi-definite. 

While  the  n2  x  n2-matrices  /  —  ^fDxADx  and  I  —  DyADy  are  sparse,  their  inverses 
are  not.  However,  the  computational  cost  of  solving  linear  systems  with  these  matrices  is 
0(n 2)  -  see  Observation 

While  the  straightforward  analogue  of  the  Peaceman-Rachford  ADI  method  in  M3  is  not 
guaranteed  to  be  stable  even  if  DXADX,  DyADy,  and  DZADZ  are  negative  semi-definite, 
stable  versions  in  M3  have  been  constructed  by  Douglas,  Rachford,  Gunn  and  others  (see, 
e.g.,  (Ml  mi  and  references  therein). 
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Observation  13.  Using  (4-14)  o,nd  the  definition  of  L-stability  in  (2.10),  it  is  clear  that 
the  Peaceman-Rachford  ADI  method  is  not  L-stable,  that  is 


lim  30)  0  0, 

lZe(z)^—oo 


(4.16) 


where  g(z)  is  the  amplification  factor.  While,  strictly  speaking,  L-stability  is  not  necessary 
for  the  solution  of  parabolic  PDEs,  strong  suppression  of  the  amplification  factor  at  infinity 
is,  that  is 

lim  g(z)  <C  1.  (4-17) 

7Ze(z)^—oo 

Therefore,  the  Peaceman-Rachford  ADI  methods  is  not  an  appropriate  choice  to  drive  the 
SDC  schemes  for  the  solution  of  parabolic  PDEs. 


However,  there  exists  a  class  of  ADI  methods  for  the  solution  of  (4.8)  based  on  the 
implicit  Euler  method  that  are  L-stable  (see,  e.g., 


While  these  methods  exhibit 
convergence  of  order  one  in  contrast  to  the  Peaceman-Rachford  ADI  methods  which  exhibit 
convergence  of  order  two,  we  accelerate  the  convergence  by  using  them  to  drive  the  SDC 
schemes. 


The  implicit  Euler-based  ADI  method  for  the  solution  of  (4.8)  is 

p*  =  <p(i)  +  hiDxADx<p* 


Its  generalization  to  M3  is 


Since  we  have 


¥>(<+1)  =  <p*  +  hiDyADy^+V. 

P*  =  +  hiDxADxip* 

V**  =tf*+  hiDyADyip** 
P{i+ 1)  =  ip**  +  hiDzADz^i+1\ 


P{i+1 )  =  {I-  hiDyADy )-\l  -  hiD.ADj-1^ 


(4.18) 


(4.19) 


m 


the  method  is  stable  provided 


/)((/-  hiDyADy)  1(I  -  hiDxADx) 


\-l 


<  1. 


(4.20) 


(4.21) 


The  stability  analysis  in  M3  is  identical. 

As  a  direct  consequence  of  the  definition  of  stability  in  (2.18)  and  L-stability  in  (2.10), 
we  have  the  following  Lemma.  See,  e.g.,  |44j  for  further  details. 


Lemma  14.  The  implicit  Euler-based  ADI  method  is  stable  if  both  DXADX  and  DyADy 
are  negative  semi- definite.  The  implicit  Euler-based  ADI  method  is  also  L-stable. 

As  a  result,  this  method  is  an  appropriate  choice  to  drive  the  SDC  schemes  for  the 
solution  of  parabolic  PDEs. 
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5  High-Order  Differentiation  Matrices 


In  order  to  develop  high  order  methods  for  solving  PDEs,  it  is  necessary  to  develop  an  ac¬ 
curate  discrete  approximation  to  the  second-derivative  operator  with  variable  coefficients, 
DXADX.  In  constructing  DXADX,  we  would  like  to  obtain  a  symmetric,  negative  definite 
matrix  for  symmetric,  negative-definite  problems  coming  from  linear,  second-order  PDEs. 
Furthermore,  such  matrix  should  permit  fast  algorithms  for  solving  associated  linear  sys¬ 
tems. 

To  construct  the  discrete  approximation  to  the  derivative  operator,  we  first  represent 
a  function  /  by  a  truncated  series  expansion  (e.g.  Fourier  or  Legendre)  and  compute  the 
derivative  by  analytically  differentiating  the  series.  The  result  is  a  linear  mapping  between 
the  vector  of  m  function  values  {/(a^)}"^1  to  the  vector  {/'(xj)}” Lq1.  Such  approach  is 
numerically  stable  provided  that  the  quadrature  nodes  are  selected  carefully.  For  example, 
if  f{x)  is  periodic,  we  use  equally-spaced  nodes.  If  f(x)  =  0  on  the  boundary  of  the  domain 
(homogenous  Dirichlet  boundary  conditions),  then  we  use  Gauss-Legendre  or,  alternatively, 
Chebyshev  nodes.  In  this  case,  the  matrix  representation  of  the  linear  mapping  between 
function  values  and  values  of  the  derivative  is  referred  to  as  the  spectral  differentiation 
matrix. 

For  details  on  the  construction  of  differentiation  matrices  for  periodic  boundary  condi¬ 
tions,  see,  e.g.,  |23l  133] .  For  details  on  the  construction  of  high-order  differentiation  matrices 
for  homogeneous  Dirichlet  boundary  conditions,  see,  e.g.,  pa  113  EH  El  S3].  We  summarize 
the  main  details  here. 


5.1  Periodic  Boundary  Conditions 

In  order  to  describe  discrete  approximations  of  the  derivative  operator  for  a  periodic  func¬ 
tion,  we  first  introduce  the  Discrete  Fourier  Transform  (DFT).  The  DFT  is  a  transformation 
T  :  Cn  — >  Cn  defined  by 

^  71—1 

fe  =  ~  Ylfie~2ni£j/n’  *  =  0,...,ra-l.  (5.1) 

U  3=0 

The  inverse  DFT  is  given  by  the  formula 

n—  1 

fj  =  ^Zfee2^j/n,  j  =  0, . . .  ,n  —  1  (5.2) 

k= 0 


(see,  e.g.,  [23,33]).  The  n  complex  numbers  {fp}  are  referred  to  as  the  Fourier  coefficients 
of  the  function  values  {/,}.  In  many  applications,  such  as  those  of  this  dissertation,  it  is 
convenient  to  view  {fj}  as  n  equally-spaced  samples  of  a  continuous,  periodic  real  function 
/  with  a  single  period  defined  on  [0,1]. 


Observation  15. 

formulae 


The  DFT  and  its  inverse  can  also  be  described  by  the  closely  related 


71—1 


h  =  lY.fje~2niij/n,  t 


3=0 


n  —  1 
2 


n  —  1 


2 


(5.3) 
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re  —  1 

fj  =  E  h^Mj/n,  j  =  0, . . .  ,n  —  1.  (5.4) 

p—  n  —  1 

2 


While  the  forms  (5.1)  and  (5.2)  are  the  standard  representation  for  the  DFT,  the  forms 
(5.3)  and  (5-4)  are  usually  preferred  in  applications  of  DFT  to  analysis  (see,  e.g.,  fTTf).  If 


xj  =  n’  J  =  0, . . . ,  77.  —  1  are  the  equally -spaced  samples  on  [0, 1],  then  (5-4)  can  be  restated 


as 


f(xj)  = 


re  — 1 
2 

E 


fee 


2ni£ti 


j  =  0,.. .  ,n  -  1. 


(5.5) 


e=- 


re  — 1 
2 


The  corresponding  trigonometric  polynomial  for  the  evaluation  of  fix)  anywhere  on  [0, 1] 
is,  therefore, 

re  — 1 

f(x)=  ±  fee2Mx.  (5.6) 

D—_  n  —  1 
2 


Obviously,  straightforward  computation  of  the  DFT  or  its  inverse  costs  0(n2)  opera¬ 
tions.  Algorithms  that  perform  such  computations  in  C(nlogn)  are  known  as  Fast  Fourier 
Transforms  (FFT)  (see,  e.g.,  [23,  33]). 

Given  a  function  f(x)  defined  as 


re  — 1 

f(x)=  ±  feeM, 

D—_  n  —  1 
2 


the  derivative  is  computed  as 


fix) 


re  — 1 
2 

E 


2t n£fee2ni£x. 


re  — 1 
2 


(5.7) 


(5.8) 


Therefore,  given  a  periodic  function  f(x),  with  a  single  period  defined  on  [0, 1],  we  compute 
its  derivative  by  first  computing  its  Fourier  coefficients  fi  via  the  FFT,  multiplying  the 
Fourier  coefficients  by  a  factor  2ttH,  and  then  applying  the  inverse  FFT  to  obtain  f'(x). 


5.2  Dirichlet  Boundary  Conditions:  Spectral  Derivative  on  a  Single  In¬ 
terval 


In  Section  3.3,  we  used  Legendre  polynomials  for  the  purpose  of  constructing  the  spectral 


integration  matrix.  We  will  now  use  the  normalized  Legendre  polynomials  {Pj}  in  (2.29) 
to  construct  the  spectral  differentiation  matrix.  Since  the  family  of  polynomials  {Pj}  is 
orthonormal,  given 

fc-i 


f(x )  =  '52<XjPj(x), 

3=0 


(5.9) 
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we  have 


a.j  =  /  f(x)Pj(x)dx. 


(5.10) 


i-i 


Lemma  16.  Suppose  that  f(x)  :  [—1, 1]  — >  M  is  given  by  the  expansion  in  (5.9)  and  f'(x)  : 


[—1,1]  — >  M  is  defined  as 

k- 1 

f  0)  =  ^2PjPj(x). 

(5.11) 

3=0 

Then 

P  =  Ba 

(5.12) 

where 

a  =  {ao,  •  •  • ,  ak- 1}  , 

P  =  {Po,  ■  ■  -,Pk~ l}T  5 

B  is  the  k  x  k-matrix  such  that 

Bij  =  ^(1)^(1)  -  -  Ktj, 

and  K  is  the  k  x  k-matrix  such  that 

Kij  =  f  Pj{x)P'j{x)dx. 


f-i 


Proof.  We  have 

After  integrating  by  parts,  we  obtain 


Pi  =  j  f'{x)Pi(x)dx. 


'-1 


Substituting  in 


and 


Pi  =  f(l)Pi(l)-f(-l)Pi(-l)-  J  f(x)Pi(x)dx 
=  /( l)Pi(l)  -  /(-l)A(-l)  -  E  ai  f1  Ppx)H(x)dx. 

3=0  " 

k- 1 

/(i)  =  E«^(1) 

3=0 
k- 1 

i=o 


the  result  follows. 


Using  (2.28),  (2.29),  and  (2.31),  a  simple  calculation  shows  that 

=  \Jl  j  +  1  y/2i  +  1 


(5.13) 

(5.14) 

(5.15) 

(5.16) 

(5.17) 


(5.18) 


(5.19) 

(5.20) 

□ 

(5.21) 


for  j  =  i  —  1,  i  —  3,  i  —  5, . . . . 
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Observation  17.  Given  quadrature  nodes  xq. 
and  denoting 

7  =  { y/wof(x0), ...,  y/wk-if(xk-i)} 


,  Xk-i,  quadrature  weights  wq,  . . . ,  Wk-i, 
T 


we  have 
where 


a  =  (ao)  ■  ■  • ,  afc-i}  , 

7  =  Va 


(5.22) 

(5.23) 

(5.24) 

(5.25) 

(5.26) 

(5.27) 

(5.28) 

If  xq,  . . .  ,Xk~ i  are  Gauss- Legendre  nodes  and  wq,  ■  ■  ■  ,Wk— l  the  corresponding  weights, 


Vij  —  y/WiPj(Xi). 

It  is  straightforward  to  see  that  the  spectral  differentiation  matrix,  D,  is  given  by 


That  is, 
where 


D  =  VBV 
9  =  Df 


-i 


g  =  {Vmf'(x0), . . . ,  y/wk-if(xk-i)}T  ■ 


then  it  is  clear  from  (5.9)  that 


k- 1 


y/w~if(xi )  =  s/uTj  ajPj(a 

3= 0 


(5.29) 


Likewise,  it  follows  from  the  discretization  of  the  integral  in  (5.10)  that 


k- 1 


Oij  =  ^  y/w~iPj(Xi)y/mif(Xi). 


i= o 


If  we  denote  by  U  the  matrix  with  elements 

U ji  —  y/Wj  Pj  (x  j  ) , 


(5.30) 


(5.31) 


then  (5.29)  and  (5.30)  imply  thatU  =  V  1 .  A  simple  calculation  shows  thatV  is  orthogonal, 
i.e.,  U  =  VT. 


The  second-order  differentiation  matrix 

D 2  =  VB2U 


(5.32) 


obtained  from  (5.26)  maps  the  values  of  the  function  f(x)  at  xq,  . . .  ,Xk-\  to  its  second 
derivative  f"(x)  at  the  same  quadrature  nodes.  In  order  to  enforce  homogeneous  Dirichlet 
boundary  conditions  within  numerical  PDE  solvers,  we  first  construct  an  augmented  second- 
order  differentiation  matrix  by  adding  quadrature  nodes  at  x  =  —  1  and  x  =  1,  i.e.,  we  use 
Gauss-Lobatto  nodes  {xq,  . . .  ,Xk-i  are  therefore  the  interior  Gauss-Lobatto  nodes).  We 
denote  by  D 2  a  particular  instance  of  the  second-order  differentiation  matrix  D 2  computed 
at  the  Gauss-Lobatto  nodes.  Since  /(— 1)  =  /( 1)  =  0,  we  denote  by  Dq  the  second-order 
differentiation  matrix  obtained  by  removing  the  first  and  last  columns  and  rows  of  D2.  Such 
construction  of  Dq  results  in  a  second-order  differentiation  matrix  that  is  symmetric  and 
negative  definite  (see,  e.g.,  mm)- 
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Observation  18.  It  is  well  known  that  differentiation  is  not  a  numerically  robust  procedure. 
Indeed,  the  maximum  eigenvalue  of  the  spectral  differentiation  matrix  D  is  of  the  order 
0(k2).  This  implies  that  the  maximum  eigenvalue  of  D2  and,  similarly,  Dq  is  of  the  order 
0(k 4).  This  presents  a  significant  problem  in  computation  as  it  limits  the  number  of  nodes 
that  can  be  used  in  the  spatial  discretization.  See,  e.g.,  m\n\ 
the  spectral  properties  of  spectral  differentiation  matrices. 


~J2  for  more  details  on 


5.3  Dirichlet  Boundary  Conditions:  Composite  Spectral  Derivative 

Since  the  spectral  differentiation  matrix  in  D  has  a  maximum  eigenvalue  of  order  0(k2) 


(see  Observation  18  above),  in  practical  computations  we  are  limited  in  the  number  of 
nodes,  k.  An  alternative  approach  is  to  fix  k  and  instead  subdivide  the  interval.  Then,  we 
represent  the  function  /  by  a  truncated  series  expansion  on  each  subinterval  and  compute 
the  derivative  by  analytically  differentiating  the  series.  It  turns  out  that  such  an  approach 
yields  a  discrete  approximation  to  the  second-derivative  operator  with  variable  coefficients, 
DXADX ,  that  is  block  five  diagonal  (or  block  tridiagonal  in  certain  implementations),  sym¬ 
metric,  and  negative  definite,  which  permits  us  to  efficiently  solve  linear  systems  with  the 
matrix  I  —  hjDxADx  (see,  e.g.,  [2]  f58j). 

Specifically,  let  [—1,1]  be  subdivided  into  2n  intervals  of  equal  size.  For  clarity,  we 
assume  that  the  number  of  points  per  subinterval  k  is  constant.  We  define  f(x)  :  [—1,1]  — >  M 
to  be 


2n  — 1  k- 1 

/(*)  =  J2 

m= 0  j= 0 


]  m  P jm  (x) 


where 


and 


Pjm  =  2  2  \  ^  1  Pj(2n(x  -  xm)  -  1) 


(5.33) 

(5.34) 

(5.35) 

for  j  =  0, . . . ,  k  —  1  and  m  =  0, . . . ,  2n  —  1,  which  has  support  on  the  interval  [xm,  xrn+\]. 
The  functions  Pjm(x)  have  disjoint  support  for  different  choices  of  rn.  A  simple  calculation 
demonstrates  that  for  a  fixed  rn,  the  functions  {Pjm(x)}  are  orthonormal  for  all  j.  This  is 
summarized  in  the  following  Lemma. 

Lemma  19.  The  family  of  polynomials  {Pjrn{x)\  is  orthonormal  on  the  interval  [—1,1]. 


xm  =  21  nm  —  1 


Hence,  given  f(x)  as  defined  in  (5.33), 


f'Km+l 


Otjm  — 


Pjm(x)f(x)dx. 


Observation  20.  Suppose  that  xq,  . . .  x^-i  are  Gauss-Legendre  nodes,  wq,. 
corresponding  weights,  and  we  denote 


(5.36) 
. ,  Wk- 1  the 


xu  =  xi  +  21  n(Xl  +  X) 


(5.37) 
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for  all  i  =  0, . . . ,  k  —  1  and  l  =  0, . . . ,  2n  —  1.  Then,  due  to  (5.33), 


2n—l k- 1 

Vuff(xu)  =  Vwi  EE  QtjmPjm  (%il )  •  (5.38) 

m=0  j= 0 


Likewise,  discretizing  (5.36)  yields 


1 


fe-i 


—  2 n  ^  ^  \/ WjPjm{Xil)y/ Wjf  (Xji) 
i= 0 

/or  all  j  =  0, ...  ,k  —  1  and  m  =  0, . . . ,  2n  —  1 . 

Introducing  the  notation 

7  i+fc.z  =  Vunf(xii) 

and 

Qj-\-k-m  =  Qjm 

where  k  is  the  number  of  points  per  subinterval,  we  observe  that 


f  =  Va 


and 


a  =  Uf 

where  V  and  U  are  the  2 nk  x  2nk  block- diagonal  matrices  with  entries 


and 


Pi+k-l,j+k-m  —  • \JlViP jm{xil ) 


U j+k-m,i-\-k-l  —  2'n  fil-1'),  Pjm  ( xil )  ■ 


(5.39) 

(5.40) 

(5.41) 

(5.42) 

(5.43) 

(5.44) 

(5.45) 


respectively.  By  construction,  we  have  U  =  V 


-i 


Given  f(x)  defined  in  (5.33),  we  would  like  to  compute  the  same  type  of  expansion  for 
f'(x).  We  start  with  the  following  lemma. 


Lemma  21.  Given  f(x)  as  defined  in  (5.33), 

2"  —  1  k- 1 


/'(7)  =  'YlfijmPjm{x), 
m= 0  ji=0 


(5.46) 


where 


k- 1 


A;  =  (/(sl+i)P<(l)  -  /(s,)Pi(-l))  -  2n  2  (5-47) 


3= 0 


and  Kij  is  defined  in  (5.16)  or,  equivalently,  in  (5.21). 
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Proof.  By  Lemma  19  we  have  that 


Pil  = 


rxi+i 


Pu(x)f'(x)dx. 


'Xl 


Integrating  (5.48)  by  parts,  we  obtain 


Pil  =  f(x)Pu(x) \l\+1  - 


rxi+ 1 


f{x)P[i{x)dx, 


(5.48) 


(5.49) 


'Xl 


and,  substituting  (5.33)  for  f{x),  arrive  at 

Pil  = 

2n-l k-l 

-EE' 

m= 0  j= 0 


-mm- 1 


<-Xl+ 1 


Pjm{x)P'u{x)dx. 


(5.50) 


I  Xl 


Due  to  (5.34),  Pjm(x)  is  only  nonzero  on  the  interval  [xi,xi+ 1]  if  m  =  l.  Therefore, 


k-l 


Pu  =  (/(s,+i)Pi(l)  -  f{xi)Pi{- 1))  -  2n  KiMi.  (5-51) 


3=0 


□ 


Observation  22.  The  value  of  f(x)  in  (5.33)  at  interior  subinterval  boundary  values  xi  for 
l  =  1, . . . ,  2”  —  2  is  not  well-defined.  In  particular,  an  interior  subinterval  boundary  value  xi 
is  shared  by  two  subintervals  with  different  Legendre  expansions.  Generally,  the  value  at  xi 
computed  using  one  Legendre  expansion  will  differ  from  the  value  computed  using  the  other 
one.  Specifically,  for  l  /  0, 


k-l 


f{xi)  =  Vzn^atjj-xPjfl) 

3=0 


from  the  left  subinterval  and 


k-l 


f{xt)  =  VzP^ajiPji-l) 

3=0 

from  the  right  subinterval.  Similarly,  for  l  2n  —  1, 

k-l 

/(®!+l)  = 


(5.52) 


(5.53) 


(5.54) 


3=0 


from  the  left  subinterval  and 


k-l 


f(xi+ 1)  =  ^^2^1+tPp-l) 

3=0 


(5.55) 


from  the  right  subinterval. 
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To  approximate  the  interior  boundary  values,  we  use  a  combination  of  the  value  from 
the  left  and  from  the  right.  In  other  words,  for  some  0  <  a,  b  <  1,  we  define 


fc-l 

f{xi+i)  =  V2 1"^  ((1  -  a)  ■  atjiPji  1)  +  a  ■  ajj+iPj(— 1)^  (5.56) 

3=0 

and 

fc-i 

f{xi)  =  ((1  -  6)  ■  +  b  ■  a3,l- lA'W)  •  (5-57) 

i=o 

The  choice  of  a  and  b  impacts  the  resulting  approximation  to  the  derivative  operator,  which 
we  discuss  further  in  Section  15.3.11 

We  summarize  the  results  in  the  following  theorem. 

Theorem  23.  Suppose  that  f(x)  is  given  by  the  expansion 

2n-l  k- 1 

f(x )  =  ^2ajmPjm(x)  (5.58) 

m= 0  j= 0 

subject  to  homogeneous  Dirichlet  boundary  conditions 


/(-!)  =  /(  1)  =  0. 


Then 


where 


Pa 


for  l  =  1, . . . ,  2"  -  2, 


2"-l fc-l 

f'(x)  =  ^2PjrnPjm(x), 

m= 0  j= 0 


fc-l 

£ 

3=0 


(5.59) 


(5.60) 


(5.61) 


Pi,  0 


fc-l 

2  ^(1  —  a)Pi(l)Pj{l)  —  Kij^J  aij,o  +  aPi(l)Pj(—l)a.j,i, 

3=0 


(5.62) 


and 


fc-i 

Pi, 2n—i  =  2"  ^  (— (1  —  b)Pi(—l)Pj(—l)  —  Kjj'j  atj, 2«-i  —  bPi(—l)Pj(l)aj^-2  (5.63) 

3=0 

with 

Dff  =  -bPP-pPpi), 

D\?  =  (1  -  a)  Pi(l)Pj{l)  -  (1  -  b)  Pp-pP^-l)  -  Kij,  (5.64) 

=  aPi(l)Pj{- 1). 

for  some  0  <  a,  b  <  1. 
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For  a  given  a  and  b  in  (5.56)  and  (5.57),  we  denote  by  B  the  2 nk  x  2 nk  matrix  such  that 


Ba  =  (3 


(5.65) 


where 

fij+k-m  =  fijrn 

and  k  is  the  number  of  points  per  subinterval.  As  a  result,  the  spectral  differentiation 
matrix  is  given  by 

D  =  VBU.  (5.66) 

Observation  24.  Given  a  sufficiently  smooth  function  f(x),  the  error  of  the  approxima¬ 
tion  of  f'(x )  constructed  in  Theorem  23  is  0(hk~l)  where  h  =  2~n  is  the  length  of  the 
subintervals  (see  m  for  details). 


5.3.1  Second-Order  Differentiation  Matrix 

The  construction  of  the  second-order  differentiation  matrix  depends  on  the  choice  of  value 
for  f(x)  at  the  interior  subinterval  boundaries.  In  our  numerical  experiments,  we  set  the 
value  at  the  interior  boundary  as  the  average  of  the  value  computed  from  the  neighboring 
intervals  (i.e.,  a  =  b  =  f  in  Theorem  23).  We  denote  the  resulting  spectral  differentiation 
matrix  Dc.  Then,  D2  is  the  second-order  differentiation  matrix  that  maps  the  values  of  the 
function  f(x)  at  quadrature  nodes  xu  to  its  second  derivative,  f"(x),  at  the  same  quadrature 
nodes. 

Our  approach  for  the  construction  of  the  second-order  differentiation  matrix  extends 
the  construction  on  a  single  interval  in  Section  5.2  and,  e.g.,  m  m  mm,  to  multiple 
intervals.  This  approach  differs  from  the  construction  used  in  [2].  Specifically,  in  order  to 
enforce  homogeneous  Dirichlet  boundary  conditions  within  a  numerical  PDE  solver,  we  first 
construct  an  augmented  second-order  differentiation  matrix  by  adding  quadrature  nodes  at 
x  =  —  1  and  x  =  1.  That  is,  on  the  left  most  boundary,  we  use  left-hand  Gauss- Radau  nodes 
and  on  the  right  most  boundary,  we  use  right-hand  Gauss-Radau  nodes.  We  compute  the 
second-order  differentiation  matrix  at  the  augmented  nodes,  which  we  denote  by  D2.  Since 
/(— 1)  =  /( 1)  =  0,  we  denote  by  D20  the  second-order  differentiation  matrix  obtained  by 
removing  the  first  and  last  columns  and  rows  of  D2. 

The  extra  nodes  are  added  to  the  outside  subintervals  in  order  to  construct  a  Legendre 
series  expansion  with  k  +  1  terms  (as  compared  to  k  terms  on  the  interior  subintervals).  On 
the  outside  subintervals,  we  project  the  Legendre  polynomials  onto  the  space  of  polynomials 
that  satisfy  the  appropriate  boundary  condition.  If  the  extra  term  in  the  Legendre  series 
expansion  is  not  added,  then  the  basis  formed  by  projecting  onto  the  space  of  polynomials 
that  satisfy  the  appropriate  boundary  condition  will  not  span  that  entire  space.  This,  in 
turn,  results  in  a  spurious  zero  eigenvalue  that  is  avoided  by  our  construction. 

Observation  25.  Our  construction  of  D20  results  in  a  second-order  differentiation  matrix 
that  is  block  five  diagonal,  symmetric,  negative  definite.  The  condition  number  of  D20  is 
0(n%.kr),  where  nx  =  2n  is  the  number  of  intervals.  While  the  analytic  condition  number 
of  the  second- order  differentiation  matrix  is  n2k2 ,  our  construction  results  in  a  condition 
number  that  is  a  factor  of  0(k2)  greater.  However,  since  k  is  0(1),  this  is  a  numerically 
useful  discretization  of  the  second-order  differentiation  operator. 
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Observation  26.  Given  a  function  f(x)  at  quadrature  nodes  {xu}li=°0’' we  apply  the 
matrix  D20  to  the  vector  f  to  obtain  approximate  values  for  f"(x)  at  the  quadrature  nodes. 
If  f(x)  is  sufficiently  smooth,  then  the  approximation  error  is  0(hk~2)  where  h  =  2~n  is 


the  length  of  the  subintervals.  This  follows  directly  from  Observation  2f 


Observation  27.  Due  to  the  block  banded  structure  of  D20,  the  computational  cost  of 
solving  linear  systems  with  the  matrix  I  —  hiD20  is  0(nxk  ■  k2),  which  is  essentially  linear 
in  the  number  of  spatial  discretization  nodes. 


5.3.2  Accuracy  of  Approximation 

The  classical  choice  of  basis  for  representing  functions  satisfying  homogeneous  Dirichlet 


boundary  conditions  on  [—1,1]  is  {sin(nf(.r  +  l))}nez-  In  Figures  5.1,  5.2,  and  5.3 


illustrate  the  accuracy  in  computing  the  second  derivative  of  the  function 


we 


7 r 


sm(b-{x  +  1)) 


(5.67) 


for  varying  b  on  the  interval  x  €  [—1,1].  Because  the  function  does  not  necessarily  satisfy 
homogeneous  Dirichlet  boundary  conditions,  we  evaluate  the  second  derivative  by  applying 
the  augmented  second-order  differentiation  matrix  D2.  In  Figures 


5.1 


5.2 


and 


5.3 


we 


construct  the  second-order  differentiation  matrix  using  64,  256,  and  1024  points,  respec¬ 
tively.  In  each  case,  we  compare  the  accuracy  of  the  second-order  differentiation  matrices 
constructed  with  varying  number  of  intervals  and  points  per  interval.  Tables  0!  and ! 
display  the  corresponding  largest  singular  values  and  the  condition  number  of  the  associated 
second-order  differentiation  matrix  with  boundary  conditions  enforced.  It  is  clear  that,  in 
each  case,  there  is  a  tradeoff  between  the  accuracy  of  the  approximation  and  the  condition 


number  of  the  matrix.  For  example,  in  modest  size  problems,  such  as  that  in  Figure  5.1 
where  64  nodes  are  used,  the  difference  between  the  condition  number  for  the  second-order 
differentiation  matrix  with  one  interval  and  with  eight  intervals  is  only  a  factor  of  about 
10.  On  the  other  hand,  the  range  of  b  for  which  the  accuracy  of  the  approximation  is  bet¬ 
ter  than  10~n  is  much  larger  for  the  second-order  differentiation  matrix  with  one  interval. 


In  larger  problems,  such  as  that  in  Figure  5.3  where  1024  nodes  are  used,  the  condition 
number  for  the  second-order  differentiation  matrix  with  one  interval  is  large.  This  results 
in  a  significant  loss  of  numerical  precision  that  necessitates  constructing  the  second-order 
differentiation  matrix  with  many  subintervals. 
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Figure  5.1:  Comparison  of  the  L2-error  for  second  derivative  of  sin(&7r^y^), 
b  €  [—35,35],  using  64-point  stencils. 


Points  Per 

Interval 

Number  of 
Intervals 

Largest 
Singular  Value 
of  Dl 

Condition 
Number  of  D2 

64 

1 

1.3107  •  lCF 

1.8923  •  105 

16 

4 

1.2749  •  105 

5.1217-  104 

8 

8 

3.9033  •  104 

1.5053  •  104 

Table  1:  Condition  number  of  second  derivative  using  64-point  stencils. 


where 
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1  .Oe+OO 


3.1e-02 


9.8e-04 


3.1e-05 


9.5e-07 


3.0e-08 


9.3e-10 


2.9e-1 1 


9.1e-13 


256  Points/ 1  Interval 
32  Points, /8  Internals 
16  Points,  ;T6  Intervals 
8  Points, ,'32  Intervals 


-150  -100  -50 


50  100  150 


Figure  5.2:  Comparison  of  the  L2-error  for  second  derivative  of  sin(67r=^y^), 
b  €  [—150, 150],  using  256-point  stencils. 


Points  Per 

Interval 

Number  of 
Intervals 

Largest 
Singular  Value 
of  Dl 

Condition 
Number  of  H2 

256 

1 

3.1264  •  10s 

4.5139-  107 

32 

8 

7.3824  •  106 

2.9846  •  10° 

16 

16 

2.0217-  106 

8.0991  •  105 

8 

32 

6.2428  •  105 

2.4056  •  105 

Table  2:  Condition  number  of  second  derivative  using  256-point  stencils. 


where 
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Figure  5.3:  Comparison  of  the  L2-error  for  second  derivative  of  sin(&7r^y^),  where 
b  €  [—650,650],  using  1024-point  stencils. 


Points  Per 

Interval 

Number  of 

Intervals 

Largest 
Singular  Value 
of  Dl 

Condition 
Number  of  D2 

1024 

1 

7.8638  •  101U 

1.1354-  1010 

32 

32 

1.1810  •  10s 

4.7741  •  107 

16 

64 

3.2348  •  107 

1.2959  •  107 

8 

128 

9.9885  •  106 

3.8490  •  10B 

Table  3:  Condition  number  of  second  derivative  using  1024-point  stencils. 


5.3.3  Why  Special  Functions? 

Although  widely  suggested,  in  many  practical  environments  the  choice  of  sine  basis  is  in¬ 
appropriate  for  the  numerical  solution  of  PDEs  with  homogeneous  Dirichlet  boundary  con¬ 
ditions.  In  particular,  the  sine  basis  forces  all  even  derivatives  of  the  function  to  satisfy 
homogeneous  Dirichlet  boundary  conditions.  For  example,  the  function 

f(x)  =  cos(7 r(x  +  1))  —  cos(27t(x  +  1))  (5.68) 

satisfies  homogeneous  Dirichlet  boundary  conditions  but  does  not  admit  an  accurate  sine 
series  representation.  Specifically,  while  the  second  derivative  of  the  sine  series  vanishes 
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on  the  boundary,  f"(x)  =  Sir2  for  x  =  —  1  and  x  =  1 .  This  means  that  the  sine  series 
representation  of  f"( x)  is  not  accurate  near  the  boundary,  resulting  in  the  so-called  Gibbs 
phenomenon.  The  sine  series  coefficients  of  f"(x)  are 


/"(x)sin(n|(x  +  1)) 


67m3  (cos  (nx)  —  1) 
(n2  —  16)  (n2  —  4) 


(5.69) 


for  odd  n.  That  is,  the  coefficients  decay  like  ~  and  are  not  absolutely  convergent.  This 
results  in  a  sine  series  representation  that  converges  (slowly)  at  every  point  except  on  the 
boundary. 

On  the  other  hand,  in  Figureph4|we  use  the  example  in  ( 5.68 )  to  investigate  the  numerical 
order  of  convergence  for  the  second-order  differentiation  matrix  (with  enforced  boundary 
conditions)  D 2  0  with  different  number  of  quadrature  nodes  per  subinterval.  While  the  error 
initially  decays  exponentially  fast,  eventually  the  error  for  each  of  the  schemes  increases 
slightly  due  to  the  growing  condition  number  of  the  second-order  differentiation  matrix. 


Figure  5.4:  Comparison  of  the  L2-error  in  computing  the  second  derivative  of  (5.68)  for 
different  order  methods  as  a  function  of  the  number  of  subintervals. 


6  Spectral  Deferred  Corrections  for  Parabolic  PDEs 


6.1  In  M1 


The  variable  coefficient  diffusion  equation  in  M1  is 


dip 

dt 


d_ 

dx 
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(6.1) 


The  Spectral  Deferred  Corrections  algorithm  for  the  solution  of  (6.1 )  proceeds  very  similarly 
to  Spectral  Deferred  Corrections  for  ODEs  (see  Algorithm  [l).  Discretizing  (6.1)  using 
spatial  discretization  nodes  xq,  . . . ,  xn-\  yields  a  system  of  OD 


is  in  (4.2), 


^  =  Dx(ADx)ip, 


(6.2) 


which  is  then  solved  using  Algorithm  [lj 


Algorithm  3  Spectral  Deferred  Corrections  for  Parabolic  PDEs  in  M1 

Compute  an  initial  approximate  solution,  ip0(xi,tj),  for  i  =  0, . . . ,  n—  1  and  j  =  0, . . . ,  m—  1 


via  the  implicit  Euler  method  (14. 


for  a  chosen  discrete  approximation  to  the  derivative 


operator  DXADX  (see  Section  5.3.1). 
For  l  =  1, ...  ,L 


1.  Compute  the  residual  function  e(xi,tj)  in  (3.2)  via  spectral  integration  (see  Section 


3.3), 


2.  Compute  5{xi,tj)  via  the  implicit  Euler  scheme 


-  see  also  Section  4.1 


3.  Update  the  approximate  solution  <pl(xi,tj)  =  ip1  1(xi,tj)  +  5(xi,tj). 


Set  the  resulting  solution  as  <p(xi,tj)  =  <pL{xi,tj). 


We  investigate  the  computational  cost  of  Algorithm  [3j  We  assume  the  number  of  subin¬ 
tervals  in  the  spatial  discretization  is  nx  and  the  number  of  points  per  subinterval  is  k ,  that 
is  n  =  nxk.  The  cost  of  constructing  the  linear  system  within  the  implicit  Euler  scheme  is 
0{nxk2)  and  the  cost  of  solving  it  is  0(nxk3)  (see  Observation  27).  In  our  implementation, 
we  use  a  version  of  LU  decomposition  designed  for  block  banded  matrices.  Since  the  LU 
decomposition  of  the  matrix  must  be  computed  at  every  time  step,  the  cost  of  constructing 
the  initial  approximation  ip°(xi,tj)  is  0(nxk3m). 

We  now  look  at  the  cost  of  a  single  run  of  spectral  deferred  corrections.  The  computation 
of  the  residual  function  e(xi,tj )  involves  applying  the  m .  x  m  spectral  integration  matrix 
to  a  vector  representing  F(t,ip(t))  at  every  point  in  space.  Hence,  the  cost  of  computing 
the  residual  function  is  0(nxkm2).  The  computation  of  5(xi,tj)  involves  solving  the  same 
linear  systems  that  appeared  in  the  construction  of  the  initial  approximation.  As  a  result, 
if  we  store  the  LU  decomposition  of  the  linear  system  at  each  point  in  time,  the  cost  of 
solving  the  system  is  0(nxk2).  This  has  to  be  done  at  every  time  step,  resulting  in  a  cost 
of  0(nxk2m).  Finally,  the  cost  of  updating  the  approximate  solution  is  0(nxkm). 

If  we  run  L  iterations  of  spectral  deferred  corrections,  then  the  cost  is  0{nxkm2L)  + 
0(nxk2mL).  Therefore,  the  overall  computational  cost  is 


0{nxk  ■  m  ■  k2)  +  0(nxk  ■  m  ■  mL )  +  0(nxk  ■  m  ■  kL ). 


(6.3) 


Since,  in  practice,  the  values  of  k,  m,  and  L  are  roughly  equal  to  the  desired  order  of  the 
scheme,  the  computational  cost  scales  linearly  in  the  number  of  discretization  nodes  in  space 


42 


( nxk ),  linearly  in  the  number  of  temporal  nodes  (m),  and  quadratically  in  the  desired  order 
of  the  scheme  ( k 2 ,  mL,  kL). 


6.2  In  M2  and  M3 

In  M2  and  M3,  the  algorithm  is  identical  to  Algorithm  [3J  except  we  replace  implicit  Euler 
method  with  the  implicit  Euler-based  ADI  method  in  (|4.18[)  and  (4.19). 

We  investigate  the  computational  cost  for  the  algorithm  in  M2 .  For  simplicity,  we  assume 
the  number  of  subintervals  in  the  spatial  discretization  is  nx  in  both  the  x-  and  y-direction 
and  the  number  of  points  per  subinterval  is  k  in  both  directions  as  well.  Thus,  the  cost 
of  constructing  the  linear  system  is  0(n2k3)  in  both  the  x-  and  y-directions.  The  cost  of 
solving  the  resulting  linear  systems  is  0(n2k^)  (see  Observation  27).  In  our  implementation 


we  used  a  version  of  LU  Decomposition  designed  for  block  banded  matrices.  Since  the  LU 
decomposition  of  the  matrix  must  be  computed  at  every  time  step,  the  cost  of  constructing 
the  initial  approximation  is  0(nxk4m). 

We  now  look  at  the  cost  of  a  single  run  of  spectral  deferred  corrections.  The  computation 
of  the  residual  function  e  involves  applying  the  mxm  spectral  matrix  to  a  vector  representing 
F(t,ip(t))  at  every  point  in  space.  Hence,  the  cost  of  computing  the  residual  function  is 
0(n2k2m2).  The  computation  of  5  involves  solving  the  same  linear  systems  that  appeared  in 
the  construction  of  the  initial  approximation.  As  a  result,  if  we  store  the  LU  decomposition 
of  the  linear  system  at  each  point  in  time,  the  cost  of  solving  the  system  is  0(nxk2)  for 
each  point  in  the  x-direction.  The  cost  is  the  same  for  each  point  in  the  y-direction.  Hence, 
the  total  cost  is  0(n2k3).  This  has  to  be  done  at  every  time  step,  resulting  in  a  cost  of 
0(n2k3m).  Finally,  the  cost  of  updating  the  approximate  solution  is  0(n2k2m). 

If  we  run  L  iterations  of  spectral  deferred  corrections,  then  the  cost  is  0(n2k2m2L)  + 
0(n2k3m,L).  Therefore,  the  overall  computational  cost  is 


0(nxkz  ■  m-  k2)  +  0(n2k 2  ■  m  ■  mL)  +  0(nxk 2  ■  m  ■  kL). 


2  ;„2 


A  similar  analysis  shows  that  the  overall  computational  cost  in  M3  is 

0(n3k3  ■  m  ■  k2)  +  0(n3k3  ■  m  ■  mL)  +  0(n3k3  ■  m  ■  kL). 


(6.4) 


(6.5) 


As  is  the  case  in  M1,  k,  m,  and  L  are  roughly  equal  to  the  desired  order  of  the  scheme.  Thus, 
the  computational  cost  scales  linearly  in  the  number  of  discretization  nodes  in  space  ( n2k 2 
in  M2  and  n3k3  in  M3),  linearly  in  the  number  of  temporal  nodes  (m),  and  quadratically  in 
the  desired  order  of  the  scheme  (k2,  mL,  kL). 


7  Numerical  Examples 

In  this  section  we  present  the  results  of  several  numerical  experiments  in  which  we  compute 
the  solution  to 

^  =  V-(a(x)V0)  (7.1) 

for  different  choices  of  boundary  conditions  and  diffusion  coefficients  a. 
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7.1  Periodic  Boundary  Conditions  with  Constant  Coefficients  in  M1 

We  compute  the  solution  for  the  diffusion  equation  with  constant  coefficients 

(7- 


with  initial  condition 

0(x,  0)  =  cos(10x)  —  0.5  sin(8x) 

and  boundary  conditions 

d  d 

0(0 ,t)  =  0(27 T,t)  and  ^0(0,  t)  =  —  (j>(2ir,t).  (7.4) 

In  this  example,  we  set  a  =  1. 

In  this  case,  the  exact  solution  is 


(7.3) 


50  520  r 

at  a,6|0'2,rl 


0(x,  t)  =  cos(10x)e  100f  —  0.5sin(8x)e 


,-64t 


(7.5) 


Figure  7.1  is  a  plot  of  the  solution  0(x,i)  in  space  for  different  points  in  time.  In  Figure 
7.2  we  plot  the  decay  of  the  L2-error  of  the  solution  at  t  =  0.1  as  the  number  of  subintervals 
in  the  temporal  domain  increases.  We  compare  schemes  of  order  4,  8,  and  12.  In  each 
case,  the  number  of  nodes  in  time  m  is  equal  to  the  order  of  the  scheme  and  the  number 
of  deferred  corrections  L  is  equal  to  m  —  1.  In  all  cases,  we  use  32  equally-spaced  spatial 
discretization  nodes. 


Figure  7.1:  Solutions  of  the  diffusion  equation  (7.2)  at  different  points  in  time. 
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Figure  7.2:  Comparison  of  the  L2-error  in  computing  the  solution  of  the  diffusion  equation 
(7.2)  for  different  order  methods  as  a  function  of  the  number  of  subintervals  in  the  temporal 
domain. 


7.2  Dirichlet  Boundary  Conditions  with  Constant  Coefficients  in 


We  compute  the  solution  for  the  diffusion  equation  with  constant  coefficients 

<9</>  d  2</>  r  i 

at  =aw  ie[0',r| 

with  initial  condition 

<j)(x,0)  =  sin(12x)  —  0.5sin(8x) 

and  boundary  conditions 

t)  =  4>(2n,  t)  =  0 

In  this  example,  we  set  a  =  1. 

In  this  case,  the  exact  solution  is  given  by 

<j>(x,t)  =  sin(12x)e~144t  —  0.5  sin(8x)e_64<. 


(7.6) 

(7.7) 

(7.8) 

(7.9) 


Figure  7.3  is  a  plot  of  the  solution  ip(x ,  t)  in  space  for  different  points  in  time.  In  Figure 


7.4  we  plot  the  decay  of  the  L2-error  of  the  solution  at  t  =  0.1  as  the  number  of  subintervals 
in  the  temporal  domain  increases.  We  compare  schemes  of  order  4,  8,  and  12.  In  each 
case,  the  number  of  nodes  in  time  m  is  equal  to  the  order  of  the  scheme  and  the  number 
of  deferred  corrections  L  is  equal  tom-1.  In  our  spatial  discretization,  we  construct  the 
second-order  differentiation  matrix  with  16  subintervals  and  16  points  per  subinterval.  From 


Figure  5.2,  this  matrix  approximates  the  exact  second  derivative  of  the  initial  condition  in 


(7.7)  with  accuracy  of  about  11  digits. 
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Figure  7.3:  Solutions  of  the  diffusion  equation  (7.6)  at  different  points  in  time. 


Figure  7.4:  Comparison  of  the  L2-error  in  computing  the  solution  of  the  diffusion  equation 
(7.6)  for  different  order  methods  as  a  function  of  the  number  of  subintervals  in  the  temporal 
domain. 
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7.3  Dirichlet  Boundary  Conditions  with  Variable  Coefficients  in  M2 

We  compute  the  solution  for  the  diffusion  equation  with  variable  coefficients 

^  =  V  •  (a(x,  y)V4>)  x,  y  £  [0,  vr]  (7.10) 

with  initial  condition 

4>(x,  y ,  0)  =  sin(10x)  (cos(2y)  —  cos(4y))  (7.11) 

and  homogeneous  Dirichlet  boundary  conditions  (i.e. ,  if  T  is  the  boundary  of  the  square 
domain  [0,7r]2,  then  ip(x,y,t)  =  0  for  all  (x,y)  £  T).  In  this  example,  we  set 


a(x,  y )  =  (1.1  +  sin(4x))(l.l  +  cos(8y)). 


(7.12) 


Figure  7.5  is  a  plot  of  the  solution  cj)(x,  y,  t)  in  space  for  a  fixed  y  =  f  and  different 


points  in  time.  Figure  7.6  is  a  plot  of  the  solution  <f>(x,  y,  t )  in  space  for  a  fixed  x  =  f  and 
different  points  in  time.  In  Figure  7.7  we  plot  the  decay  of  the  L2-error  of  the  solution 


at  t  =  0.1  as  the  number  of  subintervals  in  the  temporal  domain  increases.  We  compare 
schemes  of  order  4,  8,  and  12.  In  each  case,  the  number  of  nodes  in  time  m  is  equal  to  the 
order  of  the  scheme  and  the  number  of  deferred  corrections  L  is  equal  tom-1.  In  this  case, 
the  error  of  a  solution  computed  with  a  given  number  of  subintervals  in  time  is  constructed 
by  comparing  that  solution  to  the  solution  computed  with  twice  as  many  subintervals. 

In  our  spatial  discretization,  we  construct  the  second-order  differentiation  matrix  in  the 
.T-direction  with  16  subintervals  and  16  points  per  subinterval.  From  Figure  [oT2j  this  matrix 
approximates  the  exact  second  derivative  of  the  initial  condition  in  (|7.11[)  for  a  fixed  y  with 
accuracy  of  about  11  digits.  We  construct  the  second-order  differentiation  matrix  in  the 


y-direction  with  8  subintervals  and  16  points  per  subinterval.  From  Figure  5.4,  this  matrix 


approximates  the  exact  second  derivative  of  the  initial  condition  (7.11)  for  a  fixed  x  with 
accuracy  of  about  11  digits. 
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Figure  7.5:  Solutions  (at  y  =  5)  of  the  diffusion  equation  (7.10)  at  different  points  in  time. 


Figure  7.6:  Solutions  (at  x 


?)  of  the  diffusion  equation  (7.10)  at  different  points  in  time. 
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4th  Order 


Figure  7.7:  Comparison  of  the  L2-error  in  computing  the  solution  of  the  diffusion  equation 
( 7.10 )  for  different  order  methods  as  a  function  of  the  number  of  subintervals  in  the  temporal 
domain. 


8  Conclusions  and  Future  Work 

We  introduce  a  new  class  of  algorithms  for  the  solution  of  parabolic  partial  differential 
equations.  We  combine  Spectral  Deferred  Correction  methods  for  the  solution  of  systems 
of  ordinary  differential  equations  with  the  implicit  Euler  method  and  Alternating  Direction 
Implicit  methods  in  M2  and  M3.  Furthermore,  we  construct  high-order  accurate  representa¬ 
tions  of  the  spatial  operator.  We  extend  the  traditional  pseudospectral  schemes  by  subdi¬ 
viding  the  entire  spatial  domain,  constructing  bases  on  each  subdomain,  and  combining  the 
obtained  discretization  with  the  implicit  SDC  schemes.  As  a  result,  we  construct  schemes 
with  arbitrary  order  of  convergence  for  the  solution  of  parabolic  PDEs.  Our  schemes  have 
CPU  time  requirements  that  are  linear  in  the  number  of  spatial  discretization  nodes  and 
linear  in  the  number  of  temporal  nodes. 

Schemes  up  to  order  12  in  both  time  and  space  have  been  extensively  tested.  We 
are  currently  in  the  processes  of  implementing  parallelized  and  adaptive  versions  of  the 
approach.  An  additional  advantage  of  ADI  schemes  in  higher  dimensions  is  that  they  can 
be  trivially  parallelized,  which  leads  to  greater  computational  efficiency.  Versions  of  this 
approach  for  the  solution  of  linear  and  nonlinear  problems  of  wave  propagation  are  being 
vigorously  pursued. 
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